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Introduction 


The  hypothalamic-pituitary-adrenal  (HPA)  axis  controls  the  body’s  “fight  or  flight”  response  through  a  series 
of  endocrine  and  immune  signals  directed  at  ensuring  immediate  survival  and  later  re-establishing 
homeostasis.  Changes  in  the  tone  of  this  response  have  been  observed  in  veterans  with  Gulf  War  Illness 
(GWI).  Studies  report  abnormal  cell  proliferation,  impaired  function  and  persistent  oxidative  stress  in 
circulating  immune  cells  of  patients.  Similarly  dysregulation  of  the  HPA  axis  includes  hypersensitivity  in 
cytokine  feedback  as  well  as  suppression  of  cortisol  and  neurotransmitters  responsible  for  mediating  innate 
and  adaptive  immunity.  This  is  further  complicated  by  the  impact  on  the  HPA  axis  of  a  myriad  of  regulatory 
interactions  both  within  and  between  i)  the  immune  system  and  ii)  the  sex-hormone  axis,  the  hypothalamic- 
pituitary-gonadal  (HPG)  axis. 

We  proposed  that  severe  physical  or  psychological  insult  to  the  endocrine  and  immune  systems  can 
displace  these  from  a  normal  regulatory  equilibrium  into  a  compromised  stable  state.  This  state  is 
characterized  by  a  self-perpetuating  inflammatory  response  that  involves  regulatory  imbalance  between  the 
HPA,  HPG  and  immune  axes.  To  explore  the  validity  of  this  hypothesis  our  objective  was  to  create 
comprehensive  engineering  models  of  endocrine-immune  interaction  dynamics  in  order  to  identify  (i) 
theoretical  failure  modes  of  the  endocrine-immune  interplay  that  align  with  GWI,  and  (ii)  promising  treatment 
strategies  that  exploit  the  naturally  occurring  stable  points  of  these  systems. 

Body. 

At  the  time  of  our  last  update  we  had  initiated  work  consistent  with  Task  3  in  the  original  Statement  of  Work 
(SOW).  This  analysis  had  focused  mainly  on  the  recovery  and  re-assessment  of  the  basic  differential 
equation  model  for  the  HPA  axis  reported  in  Ben-Zvi  et  al. ,  2009  [1],  We  identified  a  number  of 
shortcomings  in  this  basic  model.  In  particular,  we  assessed  the  natural  frequencies  supported  by  the  basic 
model  and  found  the  latter  did  not  naturally  support  circadian  rhythm  without  an  external  forcing  function.  In 
addition,  analysis  of  limited  experimental  data  in  female  subjects  with  a  closely  related  illness,  chronic 
fatigue  syndrome  (CFS)  [2],  suggested  that  a  linear  model  of  the  adrenal  gland  might  be  limited  in  terms  of 
fidelity.  The  apparent  need  for  increased  dynamic  complexity  created  an  important  challenge  to  the 
continued  use  of  conventional  differential  equation  modeling,  as  significant  gaps  in  the  required  parameter 
estimates  exist  in  the  literature.  To  circumvent  this  problem  we  critically  re-assessed  our  approach.  Rather 
than  rely  on  questionable  estimates  of  dynamic  parameters  we  have  shifted  our  approach  to  one  that 
exploits  a  much  larger  body  of  information,  namely  known  physiological  and  biochemical  connectivity.  The 
approach  is  inspired  from  circuit  theory  and  is  described  below.  Importantly,  this  approach  supports  the 
seamless  integration  of  kinetic  information  wherever  available,  be  it  simple  sequential  precedence,  relative 
time  scale  or  detailed  dynamics.  This  shift  in  paradigm,  together  with  the  concerted  efforts  of  Dr.  Craddock 
who  joined  in  January  of  this  year,  has  allowed  us  to  continue  moving  forward  on  Task  3  and  make 
substantial  progress  on  Task  2,  4  and  5  of  the  SOW  with  the  added  help  of  Mr.  Paul  Fritsch  (research 
intern). 

1 .  Re-assessment  and  adaptation  of  modeling  paradigm  (Task  1).  The  core  concept  of  the  approach 
we  have  used  is  connectivity.  There  is  a  substantial  body  of  anatomical  and  biochemical  data  for  many 
biological  systems  describing  the  connectivity  between  components,  however  data  regarding  the  precise 
stoichiometry  and  kinetics  of  these  interactions  in  humans  is  extremely  limited.  Existing  models  rely  heavily 
on  animal  data  as  a  source  of  kinetic  parameters,  or  adopt  general  order  of  magnitude  estimates  when  this 
data  is  lacking.  To  circumvent  this  issue  and  draw  on  this  rich  body  of  physiological  and  biochemical  wiring 
diagrams,  we  have  adopted  the  discrete  logical  network  methodology  proposed  originally  by  Thomas  et  al. 
[3,4]  and  developed  further  by  Mendoza  and  Xenarios  [5].  Analysis  of  these  networks  makes  it  possible  to 
identify  the  number  and  type  (e.g.  oscillatory,  etc...)  of  resting  states  as  well  as  their  molecular  and  cellular 
profile  without  detailed  knowledge  of  response  dynamics.  We  have  used  this  methodology  to  construct  2 
complementary  models  at  different  scales. 


2.  An  integrated  model  of  HPA-HPG-immune  interaction  (Task  3,  Task  5).  Here  we  extended  our 
previous  analysis  of  HPA  axis  dynamics  [1]  by  including  feed-forward  and  feedback  interactions  with  sex 
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hormone  regulation  and  immune  response.  A  circuit  model  was  constructed  that  linked  state  variables 
across  the  HPA  axis  with  hypothalamic-pituitary-gonadal  (HPG)  function  in  both  men  and  women,  as  well  as 
innate  (HR)  and  adaptive  (AIR)  immune  response.  The  latter  was  described  by  the  expression  of  specific 
subsets  of  cytokines.  Details  of  this  analysis  are  described  in  Appendix  A  in  the  form  of  a  completed 
manuscript  that  is  currently  under  internal  review  by  collaborating  co-authors  [6],  In  brief,  interactions 
between  the  HPA,  HPG  and  immune  systems  were  described  as  a  circuit  model  consisting  of  1 1  state 
variable  nodes  where  each  node  can  assume  one  of  three  states  at  any  point  in  time:  -1  (inhibited),  0 
(nominal)  and  +1  (elevated).  In  this  model  the  HPA  axis  was  described  by  5  state-variable  nodes: 
corticotropin-releasing  hormone  (CRH),  adrenocorticotropic  hormone  (ACTH),  cortisol  (CORT)  and 
cytostolic  glucocorticoid  receptors  (GRs),  which  unlike  membrane  bound  receptors,  dimerize  (GRD).  As  a 
first  coarse  model,  immune  function  was  described  in  terms  of  innate  immune  response  (HR)  and  adaptive 
immune  response  (AIR).  The  state  of  activation  of  the  HR  node  was  meant  to  capture  changes  in  cytokines 
produced  by  the  innate  response  including  the  pro-inflammatory  cytokines  interleukin  (IL)  1,  6,  8  and  tumor 
necrosis  factor  alpha  (TNF-a).  Also  represented  by  the  HR  node  were  IL-12  and  IL-15,  two  mediators  of  NK 
cell  function,  as  well  as  IL-23,  an  important  inflammatory  signal  contributing  to  the  Th17  response. 

Cytokines  represented  by  the  AIR  node  include  IL-2,  IL-4,  IL-5,  IL-13,  interferon-gamma  (IFN-y),  and  TNF-p. 
As  IL-10  is  produced  by  both  the  HR  and  AIR,  it  was  represented  as  a  separate  state-variable  node.  HPG 
function  was  described  by  the  levels  of  gonadotropin-releasing  hormone  (GnRH),  of  luteinizing  homone  (LH) 
as  well  as  testosterone  (TEST)  in  males  (Figure  1A)  and  estradiol  (EST)  in  females  (Figure  IB).  The  effects 
of  gender  merit  special  attention.  While  testosterone  (TEST)  exhibits  an  inhibitory  effect  on  all  levels  of  the 
HPA  axis,  estrogen  (EST)  and  progesterone  can  stimulate  or  suppress  HPA  activity  depending  on  the 
phase  of  menstrual  cycle.  It  remains  unclear  whether  this  inhibition/stimulation  of  the  HPA  axis  occurs  in 
concert  with  the  inhibition/stimulation  of  the  HPG.  As  a  result,  these  cases  were  explored  as  separate 
alternative  models  of  the  HPA-HPG  interaction  in  female  subjects. 

In  this  formalism,  the  overall  state  of  the  system  is  described  by  a  vector  of  discrete  values,  either  -1,  0  or 
+1 ,  expressed  at  each  state-variable  node  (Eq.  1 ).  This  evolves  over  time  and  at  each  time  step  the  next 
state  value  x,(f+ 1)  for  the  /'th  node  in  the  system  is  determined  from  that  node’s  current  state  as  well  as  a  set 
of  balanced  ternary  logic  statements  based  on  the  state  values  of  the  neighboring  input  nodes  and  their 
mode  of  action  (i.e.  activate  or  inhibit).  The  logic  statements  that  determine  a  state  transition  are  expressed 
as  follows  (Eq.  2-4): 


x(t)  =  (xl(t),x2(t)...xN(t)) 


(1) 


x(f  +  l)  = 


(4  (0  v  4  (0—v  4  (0)0(4  (0  v  4  (0  -  v  4(0) 


(x*(0v4(0...v4(0) 

-(4(0  v  4(0- v  4(0) 


(2) 

(3) 

(4) 


where  the  0,  v,  and  -  symbols  are  ternary  HIGH/LOW  PASS,  OR  and  NOT  operators,  x,jA  is  the  state  of  the 
/th  node’s  /h  activator,  xik  1  is  the  state  of  the  /',h  node’s  kth  inhibitor.  Equation  (2)  is  used  when  the  node 
possesses  both  activators  and  inhibitors,  Equation  (3)  when  the  node  has  only  activators  and  Equation  (4) 
when  the  node  has  only  inhibitors.  The  number  of  state  variables  determines  the  total  number  of  system- 
wide  states  such  that  a  model  of  N  state  variables  possesses  3N  states.  As  a  result  the  number  of  total 
system-wide  states  increases  rapidly  as  new  state  variable  elements  are  added. 


These  calculations  were  encoded  into  a  rapid-prototyping  Python  script  that  was  used  to  search  the  above- 
mentioned  network  for  stable  equilibrium  states.  Analysis  time  scaled  approximately  with  N  state  variables 
as  0(N2).  Results  of  these  simulations  can  be  summarized  as  follows: 

•  In  a  first  validation  of  the  methodology,  a  discrete-network  equivalent  of  the  differential  equation 
model  used  previously  [1]  was  tested  and  the  same  two  stable  states  were  recovered.  These  were  a 
normal  resting  state  and  a  hypocortisolic  state. 
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Male  subjects.  Inclusion  of  basic  immune  function  and  sex  hormone  regulation  by  the  HPG  axis  with 
HPA  function  (HPA-GR-lmmune-HPG  model)  in  male  subjects  resulted  in  the  emergence  of  3  stable 
equilibrium  states.  Once  again  the  first  state  was  that  of  normal  health.  The  second  equilibrium  state 
was  characterized  by  low  levels  of  ACTH  and  elevated  expression  of  the  glucocorticoid  receptors  GR1 
and  GRD1.  The  third  state  corresponded  to  a  hypercortisolic  condition  with  high  glucocorticoid  receptor 
expression  (GR1,  GRD1)  and  adaptive  immune  activity  (AIR),  accompanied  by  a  suppression  of  the 
HPG  axis,  innate  immune  function  (HR)  and  IL-10  concentration. 

Female  subjects.  In  the  specific  case  of  positive  feedback  along  the  HPG  axis  and  suppressive 
interaction  with  the  HPA  axis,  the  HPA-GR-lmmune-HPG  model  for  female  subjects  (Figure  1 B) 
supported  6  steady  states.  In  addition  to  3  states  equivalent  to  those  obtained  for  the  male  subjects,  we 
found  new  steady  states  that  corresponded  to  low  CRH,  ACTH,  and  CORT,  with  nominal  immune 
function  accompanied  by  elevated  GnRH,  LH/FSH,  and  EST.  This  combination  occurred  at  each  of  the 
three  low,  nominal  and  high  values  for  GRD1/GR1  expression.  Note  that  a  stable  steady  state 
characterized  by  low  cortisol  levels  was  found  only  for  female  subjects. 

Alignment  with  experimental  data.  To  validate  these  results  the  predicted  steady  states  were  first 
compared  to  steroid  and  cytokine  levels  recorded  in  male  Gulf  War  veterans  with  GWI  and  healthy 
veterans  (HCs)  as  part  of  a  sister  study  [7]  (Figure  2A).  To  compare  against  model  predictions  steroid 
and  cytokine  levels  significantly  higher  in  GWI  were  assigned  a  value  of  +1,  significantly  lower  levels 
were  assigned  -1,  and  those  showing  no  significant  change  were  assigned  0.  Experimental  cytokine 
measurements  were  scored  and  aggregated  into  cytokine  sets.  Together  with  hormone  levels,  these 
data  supported  experimental  estmates  for  5  of  the  1 1  state  variables  modeled.  Comparison  of  these 
state  variables  with  the  predicted  steady  state  (SS2)  described  by  low  TEST,  high  CORT,  high  AIR 
values,  showed  a  match  in  3  out  of  5  variables  for  the  GWI  condition  (Figure  2A).  As  experimental 
measures  for  ACTH,  GR1  and  GRD1  were  not  available,  steady  state  SSI  and  the  nominal  steady  state 
SSO  (HC)  could  not  be  distinguished  and  validated  separately  from  one  another. 


As  a  much  greater  proportion  of  women  than  men  are  affected  by  CFS,  we  compared  the  predicted 
steady  states  identified  with  the  female  model  to  experimental  data  collected  under  two  compatible 
studies  [8,9],  Taken  together,  experimental  measurements  indicated  that  female  CFS  patients  presented 
with  high  EST,  hypocortisolism,  and  with  a  spectrum  of  altered  cytokine  levels.  We  found  no  agreement 
with  the  low  EST,  high  CORT  and  high  AIR  steady  of  steady  state  SS2,  which  closely  resembled  the 
GWI  profile.  However,  there  was  a  complete  alignment  in  all  5  state  variables  between  measured  levels 
of  endocrine-immune  markers  in  female  CFS  subjects  and  the  3  steady  states  characterized  by  high 
EST,  low  CORT  and  nominal  immune  function  (SS3,  SS4  and  SS5  in  manuscript;  Appendix  A). 


3.  A  detailed  model  of  immune  circuitry  (Task  2,  Task  5).  Based  on  the  work  of  Folcik  et  al.  (2007, 
2011)  [10,1 1]  and  an  extensive  review  of  recent  literature,  we  constructed  a  wiring  diagram  that  describes 
cytokine  signaling  between  immune  cell  populations  using  the  above-mentioned  methodology  (Figure  3) 

[12]  (see  excerpts  of  manuscript  in  preparation;  Appendix  B).  The  specific  cytokines  supporting  this  model 
included  interleukin  (IL)-1,  IL-2,  IL-4,  IL-5,  IL-6,  IL-8,  IL-10,  IL-12,  IL-13,  IL-23,  IL-27,  interferon  (IFN)-y,  and 
tumor  necrosis  factor  (TNF)-a.  In  accordance  with  Folcik  et  al.  (2011)  and  to  improve  computational 
efficiency,  cytokines  were  grouped  according  to  their  dominant  action  into  either  a  monokine  (MK)  or 
cytokine  (CK)  group.  The  MK  groups  represented  the  cytokines  released  primarily  by  monocytes  (ie:  DCs  or 
dendritic  cells)  and  the  CK  groups  represented  the  cytokines  released  by  lymphocytes  (ie  NK  cells,  Th  cells, 
and  CTLs).  MK1  and  CK1  were  designed  as  pro-inflammatory  cytokine  groups,  whereas  MK2  and  CK2 
were  specified  as  anti-inflammatory  cytokine  groups.  See  Table  4  in  appendix  B  for  the  cytokine  groups 
definitions  used.  It  should  be  noted  that  the  MK2  node  in  the  final  model  also  represents  the  anti¬ 
inflammatory  dendritic  cells  (DC2).  Since  the  DC2  node  had  only  one  input  (Infection)  and  only  one  output 
(MK2),  these  2  nodes  were  lumped  into  a  single  node  consisting  of  the  MK2  cytokines. 

In  order  to  capture  interaction  with  the  endocrine  system  and  allow  integration  into  the  HPA-HPG-immune 
model  described  in  Section  1,  this  detailed  immune  model  also  included  the  hormone  inputs  of  cortisol  and 
testosterone.  Cortisol  typically  acts  to  suppress  NK  cell  activity,  Thl  activity  and  inflammatory  dendritic  cell 
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(DC)  activity.  Conversely  it  has  been  demonstrated  that  the  pro-inflammatory  cytokines,  predominantly  IL-1, 
IL6,  and  TNF-a,  generally  activate  the  HPA  axis  at  the  hypothalamic  level  by  inducing  CRH  secretion 
leading  to  increased  release  of  ACTH  and  eventually  cortisol.  With  regard  to  the  male  HPG  axis,  studies 
have  shown  that  testosterone  and  other  androgens  exert  a  suppressive  effect  on  the  IL-6  gene. 
Testosterone  is  also  reported  to  enhance  the  Thl  response  and  activate  cytotoxic  T  lymphocytes  (CTLs). 
Cytokines  can  also  exert  negative  feedback  on  the  male  HPG  axis.  It  has  been  shown  that  receptors  for  the 
pro-inflammatory  cytokines,  IFN-y  and  TNF-a,  exist  on  male  leydig  cells  and  act  to  decrease  the  production 
of  testosterone.  As  well,  TNF-a  typically  down-regulates  the  release  of  GnRH  in  the  hypothalamus  and  LH 
in  the  pituitary  gland,  eventually  causing  a  decrease  in  testosterone  levels.  As  described  in  the  model  in 
Section  1,  testosterone,  via  its  androgen  receptor,  will  inhibit  activity  of  the  HPA  axis  while  cortisol  will 
suppress  the  HPG  axis,  inhibit  the  effects  of  testosterone,  and  down-regulate  androgen  receptor 
expression. 

These  interactions  along  with  an  updated  version  of  the  cytokine  signaling  network  proposed  by  co¬ 
investigator  Folcik-Nivar  have  been  integrated  into  an  extended  model  of  immune  signaling  complete  with 
the  effects  of  stress  and  sex  hormones  for  male  subjects  only  at  this  time  (Figure  3).  A  model  for  female 
immune-endocrine  signaling,  similar  in  level  of  detail,  is  being  constructed.  As  before,  a  first  version  of  this 
model  for  the  male  system  was  analyzed  producing  results  that  can  be  summarized  as  follows: 

•  Stable  immune  response  modes  in  male  subjects.  Application  of  this  discrete  dynamical  analysis  to 
the  detailed  endocrine-immune  network  yielded  three  different  stable  attractors.  In  addition  to  the 
normal  healthy  homeostatic  point  a  first  alternate  attractor  (SSI)  was  characterized  by  lower  than 
normal  levels  of  MK1 ,  MK23,  MK27,  Thl  and  DC1  activation,  accompanied  by  higher  than  normal  MK2, 
CK2,  and  Th2  activation  with  all  other  markers  at  nominal  levels.  This  attractor  represents  an 
upregulated  anti-inflammatory  response  since  the  anti-inflammatory  cytokines  (MK2  and  CK2)  and 
immune  cells  (Th2)  are  high,  while  all  other  markers  remain  either  low  or  nominal.  The  second  alternate 
attractor  (SS2)  was  characterized  by  low  relative  levels  of  NK  and  Thl  activity,  low  testosterone 
combined  with  elevated  CK1,  CTLs,  and  cortisol  with  all  other  markers  at  nominal  levels.  This  attractor 
represents  a  hypercortisolic  and  hypoandrogenic  inflammatory  state. 

•  Alignment  with  experimental  data.  Once  again,  these  predicted  steady  states  were  compared  with 
experimental  data  used  in  Section  1  [7-9].  The  added  detail  made  available  by  this  model  allowed  for 
comparison  of  8  markers  instead  of  5.  We  found  in  GWI  that  expression  in  7  of  the  8  clinical  marker 
sets  were  in  alignment  with  the  hypercortisolic/  hypoandrogenic  inflammatory  attractor  (SS2)  while  none 
of  these  supported  alignment  with  the  anti-inflammatory  attractor  (SSI )  (Figure  4A).  Comparing  clinical 
CFS  data  to  the  predicted  stable  attractors,  we  found  expression  in  5  out  of  8  clinical  marker  sets  were 
aligned  with  the  anti-inflammatory  attractor  (SSI),  and  6  out  of  8  clinical  marker  sets  were  in  alignment 
with  the  hypercortisolic/hypoandrogenic  inflammatory  attractor  (SS2)  (Figure  4B).  These  alignments  are 
best  visualized  with  a  Sammon  projection  of  the  hamming  distances  between  the  clinical  data  and  the 
predicted  model  stable  states  (Figure  5). 

Collectively  these  simulations  of  known  endocrine-immune  circuitry  support  the  existence  alternate 
homeostatic  regimes,  some  of  which  overlap  substantially  with  observed  immune  and  endocrine  status  in 
male  GWI  and  female  CFS  subjects.  Such  overlap  with  naturally  occurring  stable  regulatory  regimes  would 
certainly  be  consistent  with  the  long-term  persistence  of  symptoms  and  the  chronic  nature  of  these 
illnesses.  This  same  characteristic  may  also  explain  why  these  illnesses  appear  in  many  ways  resistant  to 
treatment. 

4.  Early  treatment  surveys  (Task  6,  Task  7).  Analysis  of  the  above-mentioned  regulatory  signaling 
circuits  not  only  provides  information  describing  the  stable  steady  states  available  to  the  system  but  also 
extensively  describes  the  ensemble  of  transitory  states  that  lead  unequivocally  to  one  steady  state  or 
another;  these  are  said  to  lie  within  that  steady  state’s  basin  of  attraction.  Importantly,  these  subsets  of 
transitory  states  will  lead  to  that  specific  stable  state  independently  of  an  individual’s  immune  and  endocrine 
response  kinetics.  This  guaranteed  convergence  to  a  healthy  equilibrium  makes  them  attractive  as  broadly 
applicable  treatment  destination  states.  In  the  design  of  minimally  invasive  interventions  our  basic  paradigm 
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is  therefore  to  identify  the  closest  transitory  state(s)  that  lie  within  the  basin  of  attraction  that  ensures  a 
return  to  normal  homeostasis. 

Based  on  this  paradigm  we  have  started  work  on  Tasks  6  and  7.  A  first  preliminary  analysis  of  the  multi¬ 
system  model  described  in  Craddock  et  al.,  2012  (Appendix  A)  suggests  that  the  effectiveness  of  a  single¬ 
agent  interventions  targeting  the  immune  or  endocrine  systems  separately  would  be  highly  dependent  on 
person-to-person  differences  in  response  kinetics.  A  joint  endocrine-immune  therapy  that  would  moderate 
free  cortisol  and  Thl  immune  activity  concurrently  appears  much  more  promising  if  we  are  to  deliver  a 
broadly  applicable  treatment  protocol.  We  are  currently  refining  these  models  as  well  as  the  treatment 
search  algorithm  to  incorporate  the  effects  of  timescale.  This  will  make  it  possible  to  take  advantage  of 
saddle  point  states  or  unstable  intermediate  states  that  lie  between  the  basins  of  attraction.  We  are 
currently  investigating  avenues  for  exploiting  broad  classes  of  kinetic  scales  that  might  make  it  possible  to 
reduce  the  treatment  complexity  even  further  and  tailor  these  interventions  to  patient  sub-groups. 

5.  Continuing  work.  Work  is  ongoing  on  the  refinement  of  the  treatment  design  algorithm  and  a 
working  document  describing  initial  results  is  being  assembled  as  a  draft  manuscript.  Central  to  this  is  the 
implementation  of  changes  allowing  for  us  to  extend  beyond  connectivity  and  exploit  differences  in  time 
scale.  A  second  major  area  of  ongoing  work  involves  the  design  of  a  circuit  model  describing  mechanisms 
of  neuroinflammation  and  neurotransmission  in  the  brain.  Each  component  model  is  being  designed  to 
connect  seamlessly  with  the  other  modules. 

Timeline.  As  described  in  the  previous  report  dated  September  30,  2011,  the  University  of  Alberta’s 
Research  services  Office  submitted  on  behalf  of  the  principal  investigator  a  request  for  a  one-year  extension 
of  the  project  term  due  to  administrative  delays.  This  request  was  reviewed  initially  by  Ms.  Strock  and  Dr. 
Phillips  of  the  DoD  (January  23,  2012)  and  we  were  asked  to  resubmit  this  request  at  a  later  date  (6-8 
months  before  end  of  project  term).  We  have  since  confirmed  with  Dr.  Rebecca  Fisher  that  this  continues  to 
be  the  correct  course  of  action  (ref.  email  from  Dr.  Fisher  dated  September  21,  2012).  In  accordance  with 
Dr.  Fisher’s  recommendation  we  plan  to  submit  a  formal  request  for  a  no-cost  extension  in  May  2013. 

Key  Research  Accomplishments. 

In  keeping  with  the  milestones  described  in  the  project  submission  initial  efforts  were  directed  at: 

•  We  have  completed  Task  1.  Based  on  work  reported  previously  we  revised  the  modeling  approach 
and  adopted  a  logic  circuit  paradigm  that  allows  us  to  draw  on  a  much  larger  body  of  anatomical  and 
biochemical  knowledge  while  also  providing  the  flexibility  to  seamlessly  incorporate  kinetic  data 
when  available.  This  approach  was  originally  deployed  in  a  high-level  Python  code  but  has  recently 
been  re-engineered  into  a  highly  efficient  C  code.  This  has  resulted  in  an  order  of  magnitude 
improvement  in  execution  speed  and  memory  usage. 

•  We  have  essentially  revised  and  completed  Task  2.  We  completed  a  careful  audit  of  the  original 
Basic  Immune  Simulator  (BIS)  model  and  a  thorough  review  of  the  most  recent  literature  related  to 
the  relevant  immune  signaling  mechanisms.  Based  on  this  we  updated  and  translated  the  BIS  into  a 
logic  circuit  model  and  have  analyzed  the  dynamic  stability  properties  of  this  regulatory  network. 

•  We  have  made  substantial  progress  in  the  central  component  Task  3.  In  this  regard  we  have 
produced  and  analyzed  a  circuit  model  combining  a  formal  representation  of  interactions  within  and 
between  the  HPA  axis,  the  HPG  sex  hormone  axis  and  basic  components  of  the  immune  system. 

•  Consistent  with  Task  4  we  have  conducted  an  analysis  of  multi-stability  properties  of  both  the  broad 
HPA-HPG-immune  model  and  the  detailed  BIS  immune  model.  Based  on  this  analysis  we  then 
compared  the  predicted  equilibrium  states  obtained  from  each  model  with  experimental  immune  and 
endocrine  data  from  male  and  female  GWI  and  CFS  subjects.  We  found  substantial  alignment  of 
GWI  with  stable  steady  states  involving  a  hypercortisolic/  hypoandrogenic  status. 

•  One  complete  manuscript  is  currently  undergoing  final  internal  review  and  another  is  in  the  final 
stages  of  completion  with  submission  expected  before  year’s  end. 


•  We  have  re-assessed  our  approach  to  treatment  design  (Task  6,  7)  and  have  begun  a  preliminary 
search  based  on  the  HPA-HPG-immune  model  described  above.  Initial  results  support  a  combined 
hormone-immune  therapy. 

Reportable  Outcomes. 

The  results  of  these  latest  analyses  are  being  presented  at  an  upcoming  closed  meeting  sponsored  by  the 
CDC  and  the  CFIDS  Association  of  America  and  held  at  the  Cold  Spring  Harbor  Laboratory’s  Banbury 
Centre  in  Long  Island,  NY  (Sep.  30  -Oct  3,  2012).  The  draft  manuscript  Craddock  et  al.,  2012,  enclosed  as 
Appendix  A,  is  completing  internal  review  by  co-authors  and  we  expect  submission  to  the  journal  PLoS 
Computational  Biology  before  Oct.  1 , 201 2.  Similarly  we  expect  the  working  document  in  Appendix  B, 

Fritsch  et  al.,  2012,  to  be  ready  for  submission  to  the  same  journal  by  year’s  end. 

Regarding  synergy  with  complementary  research  efforts,  these  findings  were  recently  used  in  the 
preparation  of  an  invited  GWIRP  Consortium  Award  as  planned  under  Consortium  Development  Award 
GW1 00070  (prime  institution  -  Wright  State  University). 

Conclusions. 

We  intend  to  file  in  May  2013  a  formal  request  for  a  no-cost  extension  of  the  project  term  due  to  a  delayed 
start.  We  have  carried  out  a  major  shift  in  paradigm  that  ensures  greater  flexibility  and  will  support  the 
necessary  breadth  of  the  project.  The  basic  algorithmic  framework  for  this  is  now  in  place.  Both  the  BIS 
model  of  immune  signaling  and  the  basic  model  of  the  HPA  axis  have  been  audited,  updated  and  extended 
significantly.  Both  models  have  now  been  translated  into  this  new  formalism.  Simulations  based  on  these 
models  have  shown  that  the  illness-specific  effects  of  gender  are  particularly  striking.  Work  continues  on  the 
refinement  of  the  intervention  design  component.  Initial  analyses  favor  the  deployment  of  a  joint  hormone- 
immune  intervention  over  strategies  that  target  these  systems  separately. 

Personnel  receiving  support  from  this  award: 

•  Travis  Carddock,  Ph.D.,  Senior  Research  Associate,  Broderick  Laboratory,  University  of  Alberta 

•  Paul  Fritsch,  undergraduate  Honors  Physiology  student,  Research  Intern,  Broderick  Laboratory, 
University  of  Alberta 
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Figure  1.  A  discrete  circuit  model  of  HPA-HPG-immune  interactions  in  male  subjects  (A)  and  female  subjects 
(B)  in  the  specific  case  of  positive  feedback  along  the  female  HPG  axis  and  suppressive  interaction  with  the 
HPA  axis  (model  c  in  the  manuscript)  [6].  Green  directed  edges  represent  an  up-regulation  of  the  target  by  the 
source  node  whereas  a  red  terminal  edge  represents  a  suppressive  action. 
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Figure  2.  Alignment  between  the  stable  states  (SSO,  SSI,  SS2)  predicted  by  the  circuit  model  and  the  discrete 
expression  (e.g.  elevated  (+1),  normal  (0),  or  suppressed  (-1))  of  5  markers  in  experimental  data  from  (A) 
male  GWI  subjects  and  (B)  female  CFS  subjects  [6], 
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Figure  3.  Wiring  diagram  for  a  circuit  model  of  the  immune  signaling  network  with  hormonal  inputs  from 
cortisol  and  testosterone  (purple).  The  nodes,  Stress,  Infection,  and  Male  HPG  axis  (grey)  were  necessary 
components  of  the  program  required  as  inputs,  but  they  were  kept  consistently  nominal  so  they  would  have  no 
effect  on  the  overall  network  as  per  the  logical  functions.  Green  connections  indicate  an  activating  influence 
and  red  connections  indicate  an  inhibiting  influence  [12]. 
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Figure  4.  Discretized  experimental  data  for  (A)  GWI  and  (B)  CFS  compared  to  the  two  predicted  alternate 
stable  steady  states  (SSI  [black]  and  SS2  [grey])  [12], 
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Sammon  Projection  of  Ternary  Hamming  Distances  Between  States 


Aggregate  Ternary  Score  1 


Figure  5.  Sammon  projection  of  the  hamming  distances  between  all  five  stable  system  states.  Green  circles 
represent  the  predicted  model  attractor  states  (SSO(HC),  SSI ,  and  SS2).  Red  circles  represent  the  aggregated 
clinical  data  for  GWI  and  CFS  [12], 
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Appendix  A: 


Craddock  TJA,  Miller  DB,  Fletcher  MA,  Klimas  NG,  Broderick  G.  Towards  an  Integrative  Model  of 
Complex  Stress-Mediated  Illnesses:  Gulf  War  Illness  and  Chronic  Fatigue  Syndrome.  2012,  Under 
internal  review. 
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Abstract 


As  a  key  component  in  the  body’s  stress  response  system,  the  hypothalamic-pituitary-adrenal  (HPA)  axis 
serves  to  orchestrate  changes  across  a  broad  range  of  major  biological  systems,  and  its  dysfunction  has 
been  associated  with  numerous  chronic  diseases  including  Gulf  War  Illness  (GWI)  and  chronic  fatigue 
syndrome  (CFS).  Though  tightly  coupled  with  other  components  of  endocrine  and  immune  function,  few 
models  of  HPA  function  account  for  these  interactions.  An  important  obstacle  has  been  the  scarcity  of  in 
vivo  kinetic  data.  Here  we  extend  conventional  models  of  HPA  function  by  including  feed-forward  and 
feedback  interaction  with  sex  hormone  regulation  and  immune  response.  This  is  accomplished  using  a 
discrete  logic  representation  based  solely  on  documented  physiological  and  biochemical  connectivity 
without  detailed  response  kinetics.  A  circuit  model  linked  state  variables  across  the  HPA  axis, 
hypothalamic-pituitary-gonadal  (HPG)  function,  both  in  men  and  women,  as  well  as  innate  (MR)  and 
adaptive  (AIR)  immune  response  described  in  terms  of  specific  cytokine  sets.  Results  indicated  that  unlike 
current  models,  detailed  representation  of  glucocorticoid  receptor  dimerization  is  not  required  to  produce 
multi-stability  when  immune  and  sex  hormone  regulation  are  considered.  Inclusion  of  these  effects  resulted 
in  at  least  two  alternate  homeostatic  regimes.  Experimental  data  in  male  GWI  subjects  showed  alignment  in 
3  of  5  endocrine-immune  markers  with  predictions  of  a  naturally  occurring  steady  state.  In  female  CFS 
subjects,  the  expression  of  all  5  endocrine-immune  markers  aligned  with  an  alternate  homeostatic  state 
predicted  by  the  model  based  on  regulatory  circuitry.  Indeed  we  found  that  hypercortisolism  may  be 
inextricably  linked  to  an  imbalance  in  the  immune  response  whereas  hypocortisolism  may  be  more  strongly 
related  to  sex  steroid  suppression  of  the  HPA  axis  and  promotion  of  HPG  function,  a  configuration 
seemingly  available  only  to  female  subjects.  Though  coarse,  these  models  may  nonetheless  support  the 
design  of  simple  and  robust  treatment  avenues. 
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Author  Summary 

The  hypothalamic-pituitary-adrenal  (HPA)  axis  links  the  brain,  endocrine  and  immune  systems  of  the  body, 
and  serves  as  a  major  regulator  of  the  body’s  response  to  stress  including  emotion,  digestion,  metabolism, 
sensory  perception,  etc....  Not  surprisingly,  dysregulation  of  the  HPA  axis  has  been  associated  with 
numerous  chronic  diseases,  of  which  Gulf  War  Illness  (GWI)  and  chronic  fatigue  syndrome  (CFS)  are  two 
examples.  While  HPA  function  is  inextricably  linked  with  that  of  its  regulatory  neighbors  few  models  account 
for  these  all-important  interactions,  as  many  of  the  key  parameters  needed  have  yet  to  be  determined 
experimentally.  Here  we  extend  the  common  models  of  the  HPA  axis  by  including  interaction  with  sex 
hormone  regulation  and  immune  response.  We  use  documented  physiological  and  biochemical  data  to 
construct  a  connectivity  network  describing  these  systems  for  both  men  and  women,  and  apply  a  discrete 
logic  to  represent  the  states  of  the  system  qualitatively  without  the  need  for  kinetic  parameters.  Results 
indicate  that  this  extended  system  can  accommodate  other  naturally  occurring  persistent  states  in  addition 
to  healthy  response  to  stress.  Two  of  these  alternate  resting  states  align  with  experimental  data  describing 
GWI  and  CFS  subjects. 
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1.  Introduction 


The  hypothalamic-pituitary-adrenal  (HPA)  axis,  a  key  component  in  the  body’s  stress  response  system, 
serves  to  articulate  changes  in  broad  range  of  homeostatic  regulators  as  a  function  of  environmental  cues. 
Such  cues  can  consist  of  both  physical  stressors  (injury,  disease,  thermal  exposure)  and  psycho-emotional 
stressors  (frustration,  fear,  fight  or  flight  decisions)  alike.  Instantiation  of  this  survival  program  is 
accomplished  through  controlled  modulation  of  the  neuroendocrine  and  immune  systems,  as  well  as  the 
sympathetic  nervous  systems  [1-3].  Considering  its  function  as  a  broad-reaching  integrator  of  major 
biological  systems,  it  is  no  surprise  that  numerous  chronic  diseases  have  been  associated  with  abnormal 
regulation  of  the  HPA  axis,  including  major  depression  [4,  5],  post-traumatic  stress  disorder  (PTSD)  [6-8], 
Alzheimer’s  disease  [9],  Gulf  War  illness  (GWI)  [10-12],  and  chronic  fatigue  syndrome  (CFS)  [13-15],  When 
compared  to  non-deployed  veterans,  Golier  et  al.,  (2007)  [10]  found  that  symptomatic  Gulf  War  veterans 
without  psychiatric  illness,  as  well  as  veterans  with  PTSD  alone,  showed  significantly  greater  cortisol 
suppression  to  dexamethasone  (DEX)  suggesting  markedly  enhanced  negative  feedback  along  the  HPA 
axis.  Further  study  by  these  same  investigators  indicated  that  this  might  be  due  to  a  significantly  attenuated 
ACTH  response  by  the  pituitary  in  veterans  with  GWI  without  PTSD  compared  to  non-exposed  subjects  and 
veterans  with  PTSD  alone  [11,  12],  A  similar  suppression  of  cortisol  response  to  DEX  was  found  in  CFS 
subjects  by  Van  Den  Eede  et  al.  (2007)  [13]  with  this  being  further  exacerbated  by  oestrogen  intake.  With 
regard  to  HPA  circadian  dynamics,  CFS  subjects  were  found  to  exhibit  significantly  increased  ACTH- 
adrenal  signaling  and  marginally  increased  inhibitory  feedback  when  compared  with  control  subjects  and 
CFS  subjects  comorbid  with  fibromyalgia  (FM)  [14,  15],  Conversely  the  pain-dominant  CFS-FM  subjects 
showed  significantly  blunted  cortisol  inhibitory  feedback.  While  evidence  such  as  this  implicates  abnormal 
regulation  of  HPA  function  leading  to  chronic  hypocortisolic  and  hypercortisolic  states  in  these  illnesses,  it  is 
not  clear  what  causes  this  dysregulation. 

Previously  we  investigated  the  possibility  that  some  of  these  pathological  states  may  coincide  with  naturally 
occurring  states  of  the  regulatory  system  in  question  [16],  These  backup  regimes  would  offer  alternate 
homeostatic  control  in  crisis  situations  but  at  the  cost  of  reduced  function.  The  existence  of  such  multiple 
stable  states  is  characteristic  of  systems  that  incorporate  feed-forward  and  feedback  mechanisms. 
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Feedforward  loops  in  biology  play  the  crucial  role  of  driving  rapid  acute  responses,  while  feedback  loops  will 
generally  limit  the  extent  of  a  response.  Both  will  also  drive  complex  dynamic  behavior,  including 
differentiation  and  periodicity  [17].  Under  standard  conditions  such  systems  tend  to  reside  in  a  given  stable 
attractor  state.  Small  perturbations  to  the  system  force  it  to  depart  temporarily  from  this  state,  only  to  be 
drawn  back  to  the  original  resting  state  once  the  perturbation  is  removed.  If  however,  the  perturbation  is  of 
significant  strength  and  duration,  the  system  may  abandon  its  normal  operating  regime  and  instead  be 
drawn  into  the  basin  of  attraction  of  a  new  alternate  resting  state.  Knowledge  of  the  system  dynamics  can 
allow  us  to  map  these  different  regimes  and  several  mathematical  models  of  the  HPA  exist  [18-26]. 
However  only  one  such  model  so  far  is  known  to  accommodate  multi-stability  in  the  dynamic  behavior  of  the 
HPA  axis.  It  does  so  via  the  addition  of  a  feed-forward  mechanism  involving  dimerization  of  the 
glucocorticoid  receptor  (GR)  complex  [27]  (Figure  1).  However,  this  model  and  the  majority  of  other  models 
do  not  extend  beyond  the  physiological  boundaries  of  the  HPA  axis  itself.  As  discussed  in  the  following 
sections,  HPA  activity  is  intertwined  with  the  behavior  of  the  hypothalamic-pituitary-gonadal  (HPG)  axis  and 
the  immune  system,  among  others,  and  this  interplay  cannot  be  ignored  when  considering  the  number  and 
nature  of  stationary  states  available  to  the  overarching  system. 

There  is  a  substantial  body  of  anatomical  and  biochemical  data  for  many  biological  systems  describing  the 
connectivity  between  elements,  the  presence  of  recurring  structural  motifs  and  functional  modules. 

However  data  regarding  the  precise  stoichiometry  and  kinetics  of  these  systems  in  humans  is  extremely 
limited.  Many  existing  models  rely  heavily  on  animal  data  as  a  source  of  kinetic  parameters,  or  adopt 
general  order  of  magnitude  estimates  when  this  data  is  lacking.  To  circumvent  this  issue  and  draw  on  this 
rich  body  of  physiological  and  biochemical  wiring  diagrams,  we  have  adopted  the  discrete  logical  network 
methodology  proposed  originally  by  Thomas  et  al.  [28,  29]  and  developed  further  by  Mendoza  and  Xenarios 
[30],  The  latter  makes  it  possible  to  identify  the  number  of  resting  states,  their  type  as  well  as  their 
molecular  and  cellular  description  without  detailed  knowledge  of  their  response  dynamics.  In  this  work  we 
use  this  method  to  extend  our  previous  analysis  of  the  complex  dynamical  behavior  of  the  human  HPA  axis 
by  including  its  regulatory  interactions  with  the  neighboring  HPG  axis  and  the  immune  system.  Based  on 
connectivity  information  alone,  we  show  that  multi-stability  is  easily  obtained  from  these  interacting  systems. 
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Moreover,  alternate  resting  states  are  identified  that  align  with  experimental  data  from  our  on-going  studies 
of  GWI  and  CFS. 

2  Methods 

2.1  A  Multi-systems  Model 

There  is  a  substantial  amount  of  physiological  data  describing  the  HPA,  HPG  and  immune  systems  as 
stand-alone  entities.  To  a  much  lesser  degree  there  also  exists  evidence  for  the  mutual  interactions 
between  these  systems.  The  following  sections  describe  the  experimental  evidence  used  to  infer  the 
topology  of  an  overarching  HPA-HPG-immune  interaction  network  (Figures  1-4). 

2.1.1  The  HPA  Axis 

Activation  of  the  HPA  axis  begins  at  the  paraventricular  nucleus  (PVN)  of  the  hypothalamus.  Specifically, 
afferents  transmitting  stress  related  information  in  the  brain  converge  on  the  medial  parvocellular  neurons  of 
the  PVN  inducing  the  release  of  several  peptides,  including  corticotropin-releasing  hormone  (CRH)  and 
arginine  vasopression  (AVP),  into  the  pituitary  hypophysial-portal  circulation.  CRH  and  AVP  act  in 
conjunction  on  membrane  bound  CRH-R1  receptors  in  the  anterior  pituitary  to  stimulate  adrenocorticotropic 
hormone  (ACTH)  synthesis,  and  its  rapid  release  into  peripheral  circulation.  ACTH  circulates  to  the  adrenal 
cortex  where  it  acts  on  the  membrane  bound  MC2-R  receptor  to  simulate  the  release  of  glucocorticoids 
(GCs)  (corticosterone  in  the  rat  and  cortisol  (CORT)  in  humans  and  nonhuman  primates).  To  regulate  the 
stress  response,  GCs  exert  negative  feedback  at  the  hypothalamus  and  pituitary  to  inhibit  further  synthesis 
and  release  of  CRH  and  ACTH,  respectively  [31].  This  is  the  standard  view  of  the  HPA  axis  utilized  in  the 
majority  of  models  (Figure  1  A).  However,  as  noted  by  Gupta  et  al.  [27]  circulating  glucocorticoids  act  via 
cytostolic  GC  receptors  (GRs),  which,  unlike  membrane  bound  receptors,  dimerize  (GRD)  and  translocate 
into  the  cell  nucleus  upon  activation  to  up-regulate  GR  synthesis  and  interact  with  other  relevant 
transcription  factors,  or  GC-sensitive  genes  (Figure  1  B).  As  all  nucleated  cells  possess  GRs,  GCs 
influence  practically  every  system  in  the  body.  As  described  below  major  systems  affected  by  GCs  include 
the  HPG  axis  and  immune  system. 
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2.1.2  The  HPG  Axis 


GCs  have  an  inhibitory  effect  on  the  HPG  axis,  a  key  component  that  controls  development  and  regulation 
of  the  reproductive  system,  at  all  levels  [32-36],  Activation  of  the  HPG  starts  from  brain  generated  pulsatile 
signals  that  stimulate  the  preoptic  area  of  the  hypothalamus  to  produce  gonadotropin-releasing  hormone 
(GnRH).  GnRH  is  secreted  into  the  pituitary  hypophysial  portal  bloodstream,  which  carries  it  to  the  pituitary 
gland,  where  it  activates  membrane  bound  GnRH-R  receptors,  resulting  in  the  synthesis  and  secretion  of 
luteinizing  homone  (LH)  and  follicle-stimulating  hormone  (FSH)  into  circulation.  These  gonadotropins  flow 
to  the  gonads  where  they  work  synergistically  to  promote  the  secretion  of  the  sex  steroids.  In  males,  LH 
binds  to  receptors  on  Leydig  cells  in  the  testes  to  stimulate  the  synthesis  and  secretion  of  testosterone 
(TEST).  In  females,  LH  activates  receptors  on  Theca  cells  in  the  ovaries  to  stimulate  the  release  of 
androstenedione,  which  is  aromatized  by  granulosa  cells  to  produce  estradiol  (EST),  and  progesterone 
(PROG).  TEST  negatively  feeds  back  on  the  HPG  to  inhibit  GnRH,  FSH  and  LH  secretion  and  synthesis 
[27].  This  feedback  mechanism  is  somewhat  more  complex  in  females  where,  depending  on  the  phase  of 
the  female  menstrual  cycle,  EST  and  PROG  can  exert  either  positive  or  negative  feedback  on  the 
production  and  release  of  GnRH  and  the  gonadotropins  [35,  37]. 

A  lesser-known  aspect  is  that  several  components  of  the  HPG  axis  exert  reciprocal  effects  on  the  HPA  axis 
[32,  33,  35],  Thus,  an  interactive  functional  cross  talk  exists  between  the  HPA  and  HPG  axes,  which  cannot 
be  ignored  when  investigating  HPA  axis  regulation  and  dysfunction.  Testosterone  exhibits  an  inhibitory 
effect  on  all  levels  of  the  HPA  [27]  (Figure  2  A),  whereas  EST  and  PROG  can  serve  to  stimulate  or  inhibit 
the  HPA  axis  depending  on  menstrual  cycle  phase,  or  phase  of  life  [33].  It  is  not  clear  whether  the  EST  and 
PROG  inhibition/stimulation  of  the  HPA  occurs  in  coordination  with  the  inhibition/stimulation  of  the  HPG.  As 
a  result,  these  cases  were  explored  as  separate  alternative  models  of  the  HPA-HPG  interaction  (Figure  2  B- 
D). 

2.1.3  A  Simplified  Model  of  the  Immune  System 

While  not  typically  considered  part  of  the  neuroendocrine  system,  the  immune  system  plays  a  very 
important  role  in  regulating  the  HPA  axis.  As  with  the  HPG  axis,  there  exists  a  mutual  crosstalk  between 
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these  two  systems  [38-40],  The  innate  immune  response  (MR)  is  regulated  by  cytokines  produced  primarily 
by  mononuclear  phagocytes,  such  as  macrophages,  and  dendritic  cells,  however  they  may  also  be 
produced  by  T-lymphocytes,  natural-killer  (NK)  cells,  endothelial  cells  or  mucosal  epithelial  cells.  These 
cytokines  include  the  pro-inflammatory  cytokines  interleukin  (IL)  -1,  IL-6,  IL-8  and  tumor  necrosis  factor 
alpha  (TNF-a),  and  can  also  include  IL-12,  a  primary  mediator  of  early  MR.  Also  included  is  IL-15,  which 
stimulates  proliferation  of  NK  cells  and  effector  T-lymphocytes,  as  well  as  IL-23,  an  important  inflammatory 
signal  contributing  to  the  Th17  response  against  infection.  Cytokines  that  regulate  the  adaptive  immune 
response  (AIR)  come  primarily  from  T-lymphocytes  that  have  recognized  a  cell-specific  antigen.  These 
include  IL-2,  IL-4,  IL-5,  IL-13,  interferon-gamma  (IFN-y),  and  TNF-|3.  IL-10,  which  has  important  anti¬ 
inflammatory  and  immunosuppressive  activities,  is  primarily  produced  by  macrophages,  but  is  also  produced  by 
T-lymphocytes,  dendritic  cells,  neutrophils,  eosinophils,  mast  cells,  B  cells,  NK  cells,  CD8+  T  cells,  and  all 
CD4+  T-helper  cells,  and  is  therefore  considered  separate  from  both  the  HR  and  AfR  [41], 

The  cytokines  of  the  MR  and  AIR  all  serve  to  stimulate  the  HPA  axis  at  all  levels  [38-40].  GCs,  in  turn,  serve 
to  suppress  the  synthesis  and  release  of  the  pro-inflammatory  cytokines  of  the  MR,  while  causing  a  shift 
from  the  inflammatory  response  to  the  anti-inflammatory  response  and  promotion  of  the  AIR  [38,  39,  42],  In 
the  simplified  model  used  here,  the  effect  of  the  key  regulating  cytokine  IL-10  is  produced  primarily  by  HR  with 
suppressive  effects  on  both  HR  and  AIR  (Figure  3). 

2.1.4  An  Integrative  Model  of  the  HPA-HPG-lmmune  System 

The  interactions  of  the  HPA,  HPG  and  immune  systems  outlined  above  provide  the  basis  for  a  simplified 
model  that  nonetheless  encompasses  a  broad  set  of  the  major  interactions  that  link  these  systems. 
Analyzing  the  components  of  this  model  in  a  piecewise  fashion  can  highlight  the  involvement  of  each 
individual  system  in  the  resulting  stable  attractor  states.  Piecewise  combination  of  the  above  HPG  and 
immune  systems  with  the  HPA  involves  10  individual  steps  (Figures  1-4  and  S6).  However,  as  the  action  of 
CORT  is  accompanied  by  up-regulation  of  GR  synthesis  all  connections  involving  CORT  can  include  the 
GR-GRD  feed-forward  mechanism.  As  such,  the  GR-GRD  feedback  loop  was  also  added  piecewise  to 
major  components  of  the  HPG  and  immune  systems  over  29  different  steps  (Figures  1-4  and  S1-S9).  At 
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each  step  the  resulting  HPA  axis  and  extended  system  was  analyzed  via  a  discrete  state  representation  to 
solve  for  the  stable  states. 

2.2  A  Discrete  State  Representation 

Following  the  methods  of  Thomas  et  al.  [28,  29],  and  more  recently  Mendoza  and  Xenarios  [30], 
connectivity  network  models  of  the  neuroendocrine-immune  system  are  represented  as  a  discrete  series  of 
interconnected  elements  with  three  states:  -1  (inhibited),  0  (nominal)  and  1  (activated).  The  current  state  of 
all  nodes  in  the  system  is  described  by  a  state  vector  x(t),  such  that: 

x(t)  =  (xl(t),x2(t)...xN(t)) 

where  xN(t)  is  the  state  of  the  Nth  node  of  the  system  at  time  t.  The  preferred  state  towards  which  the 
system  evolves  in  the  next  time  increment  is  described  by  the  image  vector  x(t  + 1) .  The  state  value  of  the 
image  vector  for  the  ith  node  in  the  system  is  determined  from  the  node’s  current  state  and  a  set  of  balanced 
ternary  logic  statements  based  on  the  state  values  of  the  neighboring  input  nodes  and  their  mode  of  action 
(i.e.  activate  or  inhibit).  The  logic  statements  that  determine  a  state  transition  are  expressed  as  follows  (Eq. 
2-4): 


xt(t  +  l)  = 


(4  (t)  v  xf2  (t) . . .  v  xfj  (0)0(4  (0  v  4  (t) . . .  v  4  (0) 

(4(0v4(0...v^(0) 

^(4(0  v  4(0- v  4(0) 


(2) 

(3) 

(4) 


where  the  (),  v,  and  -  symbols  are  ternary  HIGH/LOW  PASS,  OR  and  NOT  operators,  4  is  the  state  of  the 
ith  node’s  jth  activator,  x[k  is  the  state  of  the  ith  node’s  A*  inhibitor.  The  ternary  operators  given  in  Equations 


(2)  through  (4)  are  described  in  further  detail  in  Supplementary  Tables  SI  -  S3.  Equation  (2)  is  used  when 
the  node  possesses  both  activators  and  inhibitors,  Equation  (3)  when  the  node  has  only  activators  and 
Equation  (4)  when  the  node  has  only  inhibitors. 
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Applying  Equations  (2)-(4)  to  each  node  in  the  network  for  a  given  state  of  the  system,  x'"(t) ,  defines  the 
image  vector  T"'(r  +  1)  for  mth  state  of  the  system.  With  T"'(f  +  1)  defined,  the  system  may  be  updated 
asynchronously  following  the  generalized  logical  analysis  method  of  Thomas  et  al.  [28,  29],  According  to 
this  method  the  ith  element  of  the  mth  state  vector  x"'(t)  is  moved  one  step  towards  its  preferred  image 
T"'0  +  1)  (e.g.  If  xm(t)  =  -1  and  xm(t  +  l)  =  1,  then  xt(t)  is  set  to  0).  Thus,  for  each  current  state  of  the 
system  there  are  potentially  several  subsequent  states  towards  which  it  may  asynchronously  evolve. 

By  analyzing  all  possible  states  of  the  system  a  temporal  sequence  of  all  logical  states  may  be  discerned. 
To  interpret  the  results  each  state  of  the  system  can  be  represented  as  a  node  in  a  state  network.  The 
evolution  from  one  state  to  a  subsequent  state  is  then  represented  as  a  directed  edge  between  the  two 
state  nodes.  Representation  of  the  system  state  trajectory  in  this  fashion  makes  it  possible  to  draw  on  the 
concepts  and  tools  of  network  theory  for  analysis  of  the  system  dynamics.  Steady  states  are  defined  as 
those  states  for  which  the  image  vector  is  the  same  as  the  current  state  vector;  in  other  words  the  state 
node  possesses  an  out  degree  of  0. 

Twenty-nine  different  piecewise  steps  of  the  HPA  axis  models  were  analyzed  to  solve  for  stable  states 
(Table  S4).  The  number  of  state  variables  determines  the  total  number  of  system-wide  states  available  to  a 
model,  such  that  a  model  of  N  state  variables  possesses  3N  states.  Due  to  this  rule  the  number  of  total 
system-wide  states  increases  rapidly  as  new  state  variable  elements  are  added.  An  in-house  Python  script 
was  used  to  analyze  the  state  network  to  search  for  states  with  an  out  degree  of  0.  Analysis  time 
approximately  scales  with  N  state  variables  as  OfN2),  with  minor  deviations  from  this  due  to  methods  of 
parallelization,  and  other  processes  not  directly  related  to  the  analysis. 

2.3  Sample  Collection,  Processing  and  Analysis 
2.3.1  CFS  Cohort 

Endocrine  measurements  in  peripheral  blood  mononuclear  cells  (PBMC)  were  obtained  from  the  Wichita 
Clinical  dataset  [43]  for  a  group  of  39  female  CFS  subjects  and  37  Healthy  controls  (HCs)  screened  for 


confounding  medical  or  psychiatric  conditions.  Diagnostic  classification  adheres  to  the  CFS  research  case 
definition  [44],  Collection  and  processing  of  PBMCs  including  microarray  hybridization  are  found  in  [43], 
Ethics  statements  and  details  of  subject  screening,  data  preprocessing,  normalization,  outlier  detection  and 
false  discovery  correction  are  available  in  [45,  46]. 

Cytokine  profiles  were  obtained  from  a  separate  study  conducted  at  the  University  of  Miami  where  CFS 
subjects  were  diagnosed  using  the  International  Case  Definition  [47,  48],  Those  subjects  presenting  with 
additional  medical  and  psychiatric  conditions  were  excluded  from  the  study,  resulting  in  a  cohort  of  40 
female  CFS  subjects  and  a  group  of  59  healthy  female  control  subjects.  Profiling  of  cytokine  concentrations 
was  performed  in  morning  blood  plasma  samples  using  an  enzyme-linked  immuno-absorbent  assay 
(ELISA)-based  assay.  Ethics  statements  and  details  of  subject  screening,  data  preprocessing, 
normalization,  outlier  detection  and  false  discovery  correction  are  available  in  Broderick  et  al.  [49], 

2.3.2  GWI  Cohort 

Cytokine  profiles  and  endocrine  measures  were  obtained  as  part  of  a  larger  ongoing  study  of  27  GWI  and 
29  HC  subjects  recruited  from  the  Miami  Veterans  Administration  Medical  Center.  Subjects  were  male  and 
ranged  in  age  between  30  and  55.  Inclusion  criteria  was  derived  from  Fukuda  et  al.  [48],  and  consisted  in 
identifying  veterans  deployed  to  the  theater  of  operations  between  August  8,  1990  and  July  31 ,  1991 ,  with 
one  or  more  symptoms  present  after  6  months  from  at  least  2  of  the  following:  fatigue;  mood  and  cognitive 
complaints;  and  musculoskeletal  complaints.  Subjects  were  in  good  health  prior  to  1990,  and  had  no  current 
exclusionary  diagnoses  [47],  Use  of  the  Fukuda  definition  in  GWI  is  supported  by  Collins  et  al.  [50],  Control 
subjects  consisted  of  gulf  war  era  sedentary  veterans  and  were  matched  to  GWI  subjects  by  age,  body 
mass  index  (BMI)  and  ethnicity.  Ethics  statements,  methods  for  subject  assessment,  blood  analysis,  and 
data  preprocessing,  normalization,  outlier  detection  and  false  discovery  correction  are  available  in  Broderick 
et  al.  [51], 

2.3.3  Statistical  Analysis 
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All  cohort  data  was  normalized  using  a  Log2  transformation.  A  one-tailed  t-test  for  unequal  sample  sizes 
with  unequal  variance  was  used  to  determine  if  marker  levels  changed  significantly  between  subject  groups 
in  the  direction  predicted  by  the  model.  Significance  levels  were  ranked  according  to  the  null  probability  p 
with  values  less  than  0.01  considered  highly  significant,  values  ranging  from  0.01  to  0.05  considered 
significant,  values  between  0.05  and  0.10  considered  marginally  significant.  Null  probability  values 
exceeding  0.10  were  interpreted  as  non-significant. 

The  above-mentioned  experimental  data  made  it  possible  to  compare  the  expression  of  five  markers, 
namely  TEST/EST,  CORT,  MR,  IL-10,  and  AIR,  in  subjects  from  these  illness  groups  with  predictions  from 
the  model.  All  experimental  data  was  translated  from  a  continuous  scale  to  a  discrete  scale;  steroid  and 
cytokine  values  found  to  be  significantly  higher  in  patient  data  were  assigned  a  discrete  value  of  1 ,  while 
those  found  to  be  significantly  lower  were  valued  at  -1 .  Those  showing  no  significant  change  were  assigned 
a  value  of  0.  In  the  case  of  aggregate  immune  state  variables  (i.e.  MR,  IL-10,  and  AIR)  the  average  of  the 
discrete  state  scores  for  the  individual  cytokines  falling  within  each  set  was  computed  and  rounded  to  the 
nearest  integer  to  provide  an  overall  discrete  score.  Section  2.1.3  defines  the  cytokines  supporting  each 
state  variable  node.  Alignment  between  the  experimental  data  and  model  estimates  was  computed  using 
only  exact  matches.  If  the  experimental  score  and  model  prediction  for  a  given  state  variable  were  opposite 
one  another  (i.e.  +1  vs.  -1 )  or  if  one  presented  an  inactive  state  while  the  other  scored  as  active  (i.e.  0  vs. 

+1  or  -1 )  an  alignment  score  of  zero  was  assigned. 

3  Results 

3.1  Stable  States  in  the  Standard  HPA  Models 

Application  of  the  discrete  state  representation  to  the  stand-alone  HPA  model  (Figure  1  A)  generated  27 
system  states,  and  failed  to  produce  multiple  stable  states  (Table  1).  This  is  consistent  with  previous 
ordinary  differential  equation  based  models  of  the  stand-alone  HPA  axis  [16-21]  [21-26],  which  did  not  show 
system  multi-stability.  Discrete  state  representation  of  the  HPA-GR  model  (Figure  1  B)  generated  243 
system  states.  Of  these,  2  system  state  nodes  possess  no  outbound  edges  and  are  stable  attractor  steady 
states  (Table  1 ).  In  the  first  steady  state,  x(t)  =  (0,0, 0,0,0) ,  all  state  variables  assume  nominal  values.  In 


28 


the  second  steady  state,  x(t)  =  (0,-1, -1,1,1) ,  state  variables  GRD1  and  GR1  are  activated  while  ACTH  and 
CORT  are  depressed.  These  states  are  consistent  with  the  ordinary  differential  equation  models  of  the 
HPA-GR  system  used  by  Gupta  et  al.  [27]  and  Ben  Zvi  et  al.  [16]. 

3.2  The  Gender  Contribution 

Inclusion  of  the  HPG  axis  with  the  stand  along  HPA  axis  (Figure  2)  generated  729  system  states  and 
yielded  only  one  steady  state  where  all  state  variables  assume  nominal  values  (i.e.  x(t)  =  (0,0, 0,0, 0,0) ), 
except  in  the  case  of  HPG  model  c  (Table  S7).  In  model  HPA-HPGc  the  inhibition  of  the  HPA  axis  by  the 
HPG,  coupled  with  the  positive  feedback  loop  along  the  HPG  produced  an  additional  steady  state, 
x(t)  =  (-1, -1,-1, 1,1,1) ,  characterized  by  low  CRH,  ACTH  and  CORT,  and  high  GnRH,  LH/FSH  and  EST. 

When  the  HPA-GR  model  is  coupled  with  the  basic  HPG  axes  6561  system  states  are  produced.  Two 
stable  steady  states  are  generated  for  HPG  models  a,  b  and  d  (Figures  S1-S2,  Table  S7).  However,  rather 
than  reproducing  the  low  CORT  seen  in  the  HPA-GR  model,  CORT  is  held  to  its  nominal  value  while  ACTH 
is  low  and  GR1 ,  GRD1 ,  GnRH,  LH/FSH  and  EST  are  high  (i.e.  x(t)  =  (0,-1, 0,1,1, 1,1,1) ).  The  loss  of  the 
hypocortisolic  state  results  from  the  effect  of  the  sex  steroids  on  CORT.  Instead  of  CORT  receiving  only  a 
single  input  it  now  receives  multiple  inputs,  influencing  its  overall  response.  In  the  case  of  model  c  (HPA- 
GR-HPGc)  5  steady  states  were  obtained  (Figures  S2  and  Table  S7).  In  addition  to  the  two  steady  states 
obtained  with  HPA-GR-HPG  models  a,  b  and  d,  there  exist  3  new  steady  states  characterized  by  low  CRH, 
ACTH  and  CORT,  and  high  GnRH,  LH/FSH  and  EST.  The  GR1  and  GRD1  values  simultaneously  alternate 
from  low  to  nominal  to  high  to  give  the  3  new  attractor  states  (i.e.  x(t)  =  (-1,-1, — 1,-1,— 1,1, 1,1) , 

m = (-1,-1, -i,o,o, u,i) ,  m = (-1,-1, -u, u,i, i) ). 


Inclusion  of  the  GR-GRD  dimerization  feedback  loop  on  the  HPG  axis  at  LH/FSH,  and  TEST/EST  (Figures 
S3  and  S4),  produced  531441  system  states  and  further  increased  the  number  of  stable  steady  states. 
Here  the  addition  of  the  GR-GRD  loop  generates  different  profiles  for  the  four  HPG  models  (Table  S7). 
Most  system  states  present  with  low  ACTH  with  high  GRD1  and  GR1 ,  or  low  LH/FSH  with  high  GRD2  and 


29 


GR2,  or  low  TEST/EST  with  high  GRD3/4  and  GR3/4,  or  with  combinations  of  these,  while  all  other  state 
variables  assume  nominal  values.  This  set  of  system  states  is  the  result  of  sensitivity  to  CORT  introduced 
by  the  GR-GRD  feedback  loop.  The  states  of  the  HPG-GR-HPGa  model,  based  on  the  male  HPG  axis 
follow  this  regime  exclusively,  while  the  other  3  models  based  on  the  female  HPG  axis  present  with  a  variety 
of  other  states  sharing  a  hypocortisolic  characteristic.  HPA-GR-HPG-GR  models  b  and  d  show  this 
hypocortisolic  state  when  the  HPG  axis  is  suppressed  in  a  model  where  it  would  normally  promote  the  HPA 
axis.  Model  HPA-GR-HPGc  however,  only  presents  this  low  CORT  state  when  the  HPG  axis  is  promoted. 
The  multiplicity  of  system  states  that  support  a  hypocortisolic  condition  is  again  due  to  the  different 
combinations  of  activated  or  depressed  GR  and  GRD  in  the  regulatory  network. 

3.3  Immune  Effects 

Discrete  representation  of  the  simple  immune  system  joined  with  the  stand-alone  HPA  axis  (Figure  3) 
produced  729  system  states,  and  2  stable  steady  states  (Table  S6):  the  nominal  steady  state  (i.e. 
x(t )  =  (0,0, 0,0, 0,0) ),  and  a  hypercortisolic  condition  with  increased  AIR  and  depressed  HR  and  IL-10  (i.e. 
x(t)  =  (0, 0,1,1, -1,-1) ).  Incorporating  the  GRD-GR  loop  in  the  model  at  the  ACTH  state  variable  on  the  HPA 
axis  (Figure  S5  A)  generated  6561  system  states,  and  yielded  the  above  2  stable  steady  states  where  all 
immune  state  variables  were  held  at  nominal  values,  and  one  additional  state  x(t)  =  (0,-1, 0,1,1, 0,0,0) 

(Table  S6).  The  loss  of  the  hypocortisolic  state  originally  found  in  the  HPG-GR  model  is  due  to  multiple 
inputs  to  CORT,  as  also  seen  in  the  HPA-GR-HPG  models.  Adding  the  GR-GRD  feedback  loops  to  the  HR 
and  the  AIR  (Figure  S5  B)  produced  531441  system  states  and  5  stable  steady  states  (Table  S6). 

However,  the  main  characteristics  of  these  states  did  not  differ  from  the  HPA-GR-lmmune  model  as  the 
additional  states  arose  from  the  multiple  combinations  of  activated  GR5/6  and  GRD5/6. 

3.4  An  Integrative  Model 

Combining  the  HPA  axis  with  the  HPG  axis  and  immune  system  (Figure  S6)  produced  19683  system  states. 
For  HPG  models  a,  b  and  d  the  steady  state  with  all  state  variables  at  nominal  value  was  produced,  as  well 
as  a  high  CORT  steady  state  with  all  levels  of  the  HPG  held  low,  high  ADAP  and  low  HR  and  IL-10.  HPA- 
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HPG-lmmune  model  c  yielded  one  additional  state  beyond  these,  with  all  levels  of  the  HPA  held  low,  all 
levels  of  the  HPG  high,  and  all  levels  of  the  immune  system  nominal  (Table  S8).  Adding  the  HPG  axis  and 
immune  components  to  the  HPA-GR  model  (Figures  4  and  S7)  generated  177147  system  states.  HPA-GR- 
Immune-HPG  models  a,  b,  and  d  yielded  3  steady  states  (see  Table  S8).  The  first  is  the  state  with  all 
variables  at  nominal  levels,  the  second  presents  with  low  ACTH,  high  GR1  and  GRD1 ,  while  the  third 
corresponds  to  a  hypercortisolic  condition  with  high  GR1,  GRD1  and  AIR,  and  a  suppressed  HPG  axis,  MR 
and  IL-10  concentration.  HPA-GR-lmmune-HPG  model  c  presented  6  steady  states,  including  the  3 
discussed  for  the  other  models.  The  remaining  3  steady  states  showed  low  CRH,  ACTH,  and  CORT,  high 
GnRH,  LH/FSH,  and  EST,  nominal  values  for  immune  variables,  and  GRD1/GR1  state  variables  that  cycled 
between  low,  nominal  and  high. 

Previously,  the  addition  of  GR  and  GRD  to  the  LH/FSH  and  TEST/EST  variable  nodes  of  the  HPG  circuit 
(Figures  S8  and  S9)  resulted  in  different  steady  state  profiles  for  the  four  different  interaction  models.  This 
addition  was  investigated  again  with  the  addition  of  the  immune  system  component.  The  HPA-GR-lmmune- 
HPG-GR  models  produced  14348907  system  states.  Models  a,  b,  and  d  of  the  HPG  all  produced  9  steady 
states  with  the  same  profile  (Table  S8),  suggesting  that  the  immune  influence  removes  the  variability 
between  the  models.  Three  general  motifs  were  observed.  The  first  motif,  seen  in  4  of  the  steady  states, 
presents  with  a  nominal  state  for  HPA  and  immune  variables,  repressed  LH/FSH  and/or  TEST/EST  and 
increased  levels  of  the  accompanying  GR  and  GRD.  The  second,  which  presented  with  nominal  immune 
activation,  low  ACTH  and  high  GR1  and  GRD1,  as  well  as  repressed  HPG  nodes  and  increased  HPG- 
related  GR  and  GRD,  accounted  for  another  4  steady  states.  The  final  steady  state  supported  high  CORT, 
GR1  and  GRD1,  low  HPG  along  with  high  GR2/3/4  and  GRD2/3/4,  high  AIR,  and  low  IL-10  and  MR.  Model  c 
produced  these  same  9  steady  states  plus  an  additional  12  conditions  all  sharing  a  single  motif  of  low  CRH, 
ACTH  and  CORT,  high  GnRH,  LH/FSH  and  EST,  and  nominal  immune  activation  (Table  S8).  These  12 
steady  states  are  the  result  of  GR1  and  GRD1  again  cycling  between  low,  nominal  and  high  values,  and  the 
HPG  axis  GR  and  GRD  variables  cycling  between  low  and  nominal  value. 


3.5  Clinical  Indicators  of  the  GWI  Steady  State 
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To  validate  these  results  the  steroid  and  cytokine  levels  recorded  in  male  Gulf  War  veterans  with  GWI  and 
HCs  were  compared  to  the  predicted  steady  states.  At  rest  CORT  levels  in  GWI  veterans  were  found 
higher  than  normal  with  marginal  significance  (p  =  0.09),  whereas  TEST  levels  at  rest  were  found  to  be 
significantly  lower  (p  =  0.03).  The  cytokines  of  the  HR  in  GWI  patients  were  significantly  decreased  for  IL-la 
and  IL-lb  (p  =  0.05  and  0.04,  respectively),  and  significantly  increased  for  IL-6  and  IL-8  (p  =  0.03  and  0.01, 
respectively),  while  for  the  remaining  majority  of  cytokines  (IL-1 2p70,  IL-15,  IL-23  and  TNF-a)  there  was  no 
significant  change  (p  >  0.10).  There  was  no  significant  difference  in  the  expression  of  IL-1 0,  used  here  as  a 
mediating  cytokine  between  the  HR  and  AIR  state  variable  nodes  (p  =  0.36).  For  the  cytokines  belonging  to 
the  AIR  aggregate  set,  GWI  veterans  showed  a  marginally  significant  increase  in  IL-2  (p  =  0.10),  a 
significant  increase  in  IFN-y  (p  =  0.02),  and  a  highly  significant  increase  in  IL-1 3  (p  <  0.01).  IL-4,  IL-5  and 
TNF-p  did  not  show  any  significant  changes  between  groups  (p  >  0.10). 

Overall,  the  experimental  results  used  here  indicate  that  male  GWI  subjects  present  with  low  TEST, 
hypercortisolism,  and  an  increased  expression  in  the  majority  of  cytokines  related  to  the  AIR.  To  compare 
against  model  predictions  steroid  and  cytokine  values  found  to  be  significantly  higher  in  GWI  patients  were 
assigned  a  discrete  value  of  1 ,  while  those  found  to  be  significantly  lower  were  valued  at  -1 .  Those  showing 
no  significant  change  were  valued  at  0.  Experimental  cytokine  measurements  were  scored  and  aggregated 
into  cytokine  sets  as  described  above.  Comparison  of  these  state  variables  showed  a  40%  match  between 
GWI  experimental  data  and  the  nominal  SSO(HC)  stable  state  predicted  by  the  HPA-lmmune-HPGa 
iteration  of  the  model.  In  the  case  the  predicted  steady  state  SS2,  described  by  low  TEST,  high  CORT,  high 
AIR  values,  there  was  a  60%  match  with  experimental  data  for  the  GWI  condition  (Figure  5).  While,  the 
model  shows  a  difference  between  the  low  ACTH,  high  GR1  and  GRD1  steady  state  SSI  and  the  nominal 
steady  state  SSO(HC),  the  experimental  markers  do  not  allow  for  these  states  to  be  distinguished.  Thus, 
measured  endocrine-immune  profiles  for  GWI  show  a  40%  match  with  predicted  steady  state  SSI. 

3.6  Clinical  Indicators  of  the  CFS  Attractor 

Comparison  of  endocrine  marker  expression  for  female  CFS  patients  and  HCs  showed  significantly  lower 
levels  of  CORT  in  blood  (p  =  0.05),  and  a  marginally  significant  increase  in  EST  levels  (p  =  0.09)  for  the 
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fatigued  group.  Among  cytokines  belonging  to  the  MR  set,  levels  of  IL-8  were  significantly  lower  in  CFS 
subjects  (p  =  0.02),  lower  with  high  significance  for  IL-1 5  (p  <  0.01 ),  significantly  higher  for  IL-1  a  and  IL-1  b 
(p  =  0.04,  0.02),  and  higher  with  high  significance  for  IL-6  and  IL-1 2p70  (p  <  0.01 ).  IL-23  and  mediating 
cytokine  IL-1 0  showed  no  significant  change.  Of  the  AIR  cytokines,  IL-4,  IL-5  and  TNF-|3  all  showed  highly 
significant  elevations  in  expression  (p  <  0.01),  whereas  IL-1 3  levels  were  significantly  lower  in  CFS  (p  = 
0.01).  IFN-y  and  IL-2  showed  no  significant  changes. 

Taken  together,  experimental  measurements  indicated  that  female  CFS  patients  presented  with  high  EST, 
hypocortisolism,  and  with  a  spectrum  of  altered  cytokine  levels.  As  in  the  case  of  GWI,  the  average 
concentration  of  all  cytokines  falling  within  state  variable  set  was  rounded  to  the  nearest  integer  and 
projected  onto  a  discrete  scale.  Comparison  showed  a  60%  match  between  CFS  experimental  data  and  the 
nominal  (SSO(HC))  stable  state  predicted  by  the  HPA-lmmune-HPGc  iteration  of  the  model.  There  is  0% 
agreement  with  the  low  EST,  high  CORT  and  high  AIR  steady  state  SS2,  which  was  found  to  closely 
resemble  the  GWI  profile.  However,  there  was  a  100%  match  between  measured  levels  of  endocrine- 
immune  markers  in  CFS  subjects  and  the  high  EST,  low  CORT,  nominal  immune  steady  states  SS3,  SS4 
and  SS5.  As  there  are  no  clinical  markers  to  separate  the  SS3,  SS4  and  SS5  model  states  they  are 
referred  to  collectively  (Figure  6).  As  with  GWI,  the  low  ACTH,  high  GR1  and  GRD1  steady  state  SSI  and 
the  nominal  steady  state  SSO(HC),  cannot  be  distinguished  using  the  experimental  data,  and  are  also 
considered  collectively.  Thus,  the  measured  CFS  endocrine-immune  profile  showed  a  60%  match  with 
predicted  steady  state  SSI . 

Discussion 

Multiple  stable  states  are  prime  characteristics  of  systems  displaying  feedforward  and  feedback 
mechanisms,  and  play  a  critical  part  in  guiding  the  complex  dynamics  observed  in  biology.  Modeling  these 
complex  systems  can  be  quite  challenging  in  the  face  of  missing  information  regarding  the  stoichiometry 
and  in  vivo  kinetics  in  humans,  however  there  is  an  abundance  of  connectivity  data  describing  the  “wiring” 
of  these  biological  circuits.  To  make  use  of  this  wealth  of  information  we  have  applied  a  discrete  state 
representation  to  the  neuroendocrine  immune  system  based  solely  on  biological  connectivity  found  in  the 
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literature  and  a  set  of  ternary  logical  rules.  When  this  approach  was  applied  to  the  conventional  stand-alone 
HPA  axis  model  the  discrete  state  representation  failed  to  produce  multistability.  This  is  consistent  both 
with  previous  ordinary  differential  equation  models  [21-26]  and  the  analysis  of  Thomas  [29],  which 
demonstrate  that  the  presence  of  positive  feedforward/feedback  loops  generates  multiple  stable  states. 
Indeed,  addition  of  the  positive  GR-GRD  dimerization  feedback  loop  generated  a  second  stable  state 
characterized  by  low  ACTH,  CORT  and  high  GR-GRD.  This  is  also  consistent  with  the  ordinary  differential 
equation  based  models  used  by  Gupta  et  al.  [27]  and  Ben  Zvi  et  al.  [16],  Interestingly,  when  the  GR-GRD 
feedback  loop  is  added  to  the  HPA  axis  in  the  HPA-HPG  models,  the  original  low  CORT  state  found  by 
Gupta  et  al.  [27]  disappears.  This  also  occurred  with  the  HPA-lmmune  model.  This  is  a  direct  result  of  the 
state  transition  logic  applied  to  the  system.  When  multiple  inputs  determine  the  level  of  CORT  all  inputs 
must  be  low  for  CORT  to  drop,  not  just  the  level  of  ACTH  as  prescribed  by  Gupta  et  al.  [27]  and  Ben  Zvi  et 
al.  [16].  The  hypocortisolic  condition  can  be  regained  when  neighboring  regulatory  systems  are  included, 
such  as  the  female  HPG  axis.  Thus,  in  simple  models  the  GR-GRD  loop  is  a  necessary  component  to 
produce  mutlistability,  however  in  more  complex  regulatory  circuits  this  feedback  loop  becomes  redundant 
and  is  not  required  to  produce  multiple  attractor  states. 

Gender  appears  to  play  an  important  role  in  supporting  a  hypocortisolic  condition.  Due  to  the  suppressive 
actions  of  the  male  gonadal  system  in  regulating  itself  and  the  HPA  axis,  a  low  CORT  state  is  never 
available  to  the  male,  while  the  ability  of  the  female  sex  hormone  system  to  be  both  an  inhibitor  and/or 
activator  readily  supports  the  presence  of  a  hypocortisolic  condition.  In  the  absence  of  the  GR-GRD 
positive  feedback  loops  on  both  the  HPA  and  HPG  axes,  the  female  gonadal  system  is  capable  of 
generating  a  low  CORT  stable  state,  when  the  HPG  axis  inhibits  the  HPA  axis  while  promoting  itself.  While 
literature  suggests  that  this  is  a  configuration  that  is  available  to  the  female  gonadal  system  during  a  certain 
phase  of  the  menstrual  cycle,  it  remains  to  be  confirmed  experimentally.  When  the  GR-GRD  feedback  loop 
is  added  to  both  the  HPA  and  HPG  axis  multiple  hypocortisolic  states  become  available  to  the  female  in  all 
wiring  configurations  save  the  case  when  the  wiring  is  as  shown  for  the  male  (i.e.  the  HPG  inhibits  both 
itself  and  the  HPA).  Interaction  with  the  immune  system  also  appears  to  play  a  significant  role  in 
determining  abnormal  CORT  levels.  While  our  models  are  very  coarse  grained  in  regard  to  the  immune 
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system,  they  do  suggest  an  important  role  for  immune  system  feedback  in  sustaining  hypercortisolic 
conditions.  In  our  models  CORT  exerts  a  suppressive  action  on  the  MR  system  while  promoting  the  AIR.  In 
turn,  positive  feedback  by  certain  components  of  the  immune  system  further  elevates  CORT  levels  leading 
to  a  hypercortisolic  steady  state.  While,  inclusion  of  the  GR-GRD  feedback  loops  on  the  HPA  and  immune 
system  did  yield  additional  steady  states,  it  did  not  result  in  any  significant  changes  to  the  profile  in  regards 
to  CORT  levels.  The  integrative  model  including  the  HPA  and  HPG  axes  and  the  simple  immune  system 
generates  a  hypercortisolic  state  in  all  4  models  of  the  HPG  axis,  and  a  hypocortisolic  state  only  in  model  c 
of  the  HPG  axis.  The  hypocortisolic  states  found  in  the  HPA-GR-HPG-GR  models  are  removed  due  to  the 
presence  of  the  stimulatory  effects  of  the  immune  system  on  the  HPA  in  all  cases  save  for  model  c  of  the 
HPG.  Again,  the  inclusion  of  the  GR-GRD  feedback  loops  on  the  HPA  and  HPG  did  yield  additional  steady 
states,  but  did  not  result  in  any  significant  changes  to  the  system  profiles. 

These  findings  suggest  that  hypercortisolism  is  inextricably  linked  to  an  imbalance  in  the  immune  response. 
This  is  consistent  with  notions  that  GWI  may  be  due  to  a  systemic  imbalance  in  immune  signaling  [51-53], 
For  example,  Skowera  et  al.  measured  intracellular  production  of  cytokines  in  peripheral  blood  and  found 
ongoing  Thl  -type  immune  activation  in  symptomatic  Gulf  War  Veterans  compared  to  healthy  counterparts 
[52],  More  recent  work  confirms  this  finding  while  also  suggesting  that  this  may  occur  in  the  more  complex 
context  of  a  mixed  Th1:Th2  response  [51].  Additionally,  Whistler  et  al.  found  evidence  of  innate  immune 
involvement  with  GWI  subjects  exhibiting  impaired  immune  function  as  characterized  by  decreased  NK 
cytotoxicity,  altered  gene  expression  associated  with  NK  cell  function,  altered  pro-inflammatory  cytokines, 
and  T-cell  ratios  compared  to  control  subjects  [53].  This  study  also  found  dysregulated  mediators  of  the 
stress  response  (including  salivary  cortisol)  further  suggesting  a  connection  between  altered  cortisol  levels 
and  immune  imbalance  in  GWI  subjects.  Conversely  conditions  involving  hypocortisolism  may  be  more 
strongly  related  to  sex  steroid  suppression  of  the  HPA  axis  and  promotion  of  HPG  function,  a  configuration 
seemingly  available  only  to  female  subjects.  This  would  suggest  that  the  hypocortisolism  seen  in  diseases, 
such  as  CFS  [54-56],  may  be  a  result  of  the  complexity  afforded  to  the  interaction  between  the  HPA  and 
HPG  in  female  subjects,  and  may  be  an  explanation  for  the  reported  prevalence  of  such  diseases  in  women 
[57-62],  Indeed,  these  authors  report  that  approximately  70%  of  observed  CFS  patients  are  women. 
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While  certainly  more  comprehensive  than  their  predecessors,  these  models  remain  relatively  coarse 
representations  of  the  interplay  between  the  endocrine  and  immune  systems.  Misalignment  with 
experimental  data,  in  particular  in  the  case  of  GWI,  may  stem  from  several  limitations.  First  the  discrete 
logic  of  the  model  is  based  on  3  states  only  making  it  impossible  to  distinguish  between  high  and  very  high 
or  low  and  very  low.  Misalignment  may  also  be  the  result  of  missing  interactions  in  the  model  regulatory 
network.  The  missing  interactions  may  be  located  within  the  immune  network  itself,  as  immune  cells  have  a 
complex  crosstalk  involving  multiple  cytokines  [63],  However,  the  missing  interactions  may  also  involve 
systems  that  are  not  currently  modeled,  such  as  the  role  of  key  neurotransmitters  linking  the  brain  and 
central  nervous  system  with  the  HPA  axis  and  the  immune  system.  For  example,  norepinephrine  and 
epinephrine  stimulate  the  p2-adrenoreceptor-cAMP-protein  kinase  A  pathway  inhibiting  the  production  of 
Thl/proinflammatory  cytokines  and  stimulating  the  production  of  Th2/anti-inflammatory  cytokines  causing  a 
selective  shift  from  cellular  to  humoral  immunity  [64,65],  Additionally,  lymphocytes  express  most  of  the 
cholinergic  components  found  in  the  nervous  system  and  may  be  stimulated  by  or  release  acetylcholine 
thus  constituting  a  cholinergic  system,  separate  from  cholinergic  nerves,  that  regulate  immune  function  [66]. 
Similarly,  as  hypocortisolic  diseases  do  not  only  occur  in  females  the  inability  of  our  model  to  predict 
hypocortisolic  states  in  males  is  an  indication  that  low  CORT  states  in  males  may  be  due  to  finer  scale 
interactions  between  the  HPA  and  HPG  axes,  or  connections  to  other  systems,  such  as  the  brain,  that  are 
not  modeled  here.  As  an  example,  neuropeptide  Y,  a  molecule  found  in  the  nervous  system  that  activates 
and  stimulates  the  stress  response,  is  not  included  in  our  model,  but  recently  has  been  shown  to  play  a  role 
in  CFS  [67], 

The  simple  models  presented  here  illustrate  the  importance  of  an  integrative  approach  to  understanding 
complex  illnesses.  Further  refinement  of  the  model  to  include  more  detailed  description  of  interactions  within 
and  between  the  HPA,  HPG  and  immune  systems  could  extend  its  applicability  to  other  illnesses  as  would 
the  incorporation  of  other  key  systems  such  as  the  brain  and  central  nervous  systems.  Yet,  even  with  the 
coarse-grained  co-regulation  networks  investigated  we  found  numerous  stable  resting  states  that  differ 
significantly  from  normal  and  were  indicative  of  complex  and  persistent  regulatory  imbalances.  Findings 
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such  as  this  support  the  use  of  an  alternate  model  for  disease,  one  which  is  not  necessarily  associated  with 
failure  of  individual  components,  but  rather  with  a  shifts  in  their  coordinated  actions  away  from  normal 
regulatory  behavior.  Response  to  exercise  and  other  stressors  has  the  potential  to  be  very  different  in  these 
new  regulatory  regimes.  This  is  something  that  we  have  observed  firsthand  in  our  work  with  human  GWI 
and  CFS  subjects  [68],  As  these  models  are  based  on  currently  documented  knowledge  of  human 
physiology  and  regulatory  biochemistry  they  are  necessarily  incomplete.  They  may  nonetheless  provide 
enough  detail  to  identify  simple  and  robust  treatment  avenues  that  would  enable  the  body  to  recall  its 
normal  regulatory  mode  and  naturally  draw  itself  back  into  a  more  normal  homeostasis. 
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Figure  1:  The  stand-alone  HPA  models.  (A)  HPA  model.  (B)  HPA-GR  model. 


Figure  2:  The  HPA-HPG  models.  (A)  HPA-HPGa.  (B)  HPA-HPGb.  (C)  HPA-HPGc.  (D)  HPA-HPGd. 

Figure  3:  The  HPA-lmmune  model. 

Figure  4:  The  HPA-GR-lmmune-HPG-GR  models.  (A)  HPA-GR-lmmune-HPGa-GR.  (B)  (B)  HPA-GR- 
Immune-HPGc-GR.  Note:  HPA-GR-lmmune-HPGb-GR  and  HPA-GR-lmmune-HPGd-GR  may  be  found 
supplementary  information  Figure  S7. 

Figure  5:  Clinical  indicators  of  the  GWI  attractor.  Comparison  of  HPA-GR-lmmune-HPGa  model  steady 
states  with  male  GWI  clinical  data  expressed  in  terms  of  model  nodes. 

Figure  6:  Clinical  indicators  of  the  CFS  attractor.  Comparison  of  HPA-GR-lmmune-HPGc  model  steady 
states  with  female  CFS  clinical  data  expressed  in  terms  of  model  nodes. 


Table  1:  Steady  States  of  Key  HPA-axis  models.  White  -  nominal  state  (0);  Green  -  high  state  (1);  Red 
low  state  (-1 );  Grey  -  N/A  to  the  model) 


Figure  1 


Figure  2 
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Figure  3 


Figure  4 
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Figure  5 
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Appendix  B: 

Excerpts  from  Methods  section,  Tables  and  Figures  from: 

Fritsch  P,  Craddock  TJA,  Smylie  AL,  Folcik  Nivar  VA,  Fletcher  MA,  Klimas  NG,  de  Vries  G,  Broderick 
G.  A  Study  of  Multiple  Homeostatic  Regimes  in  a  Discrete  Logic  Model  of  Immune  Signaling.  2012.  In 
preparation. 
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Excerpt  from  Methods 


A  Theoretical  Immune  Signaling  Network 

Based  on  the  work  of  Folcik  et  al.  (201 1 ),  a  wiring  diagram  was  constructed  for  the  interactions  between 
immune  cells  and  their  cytokines  [Figure  1],  Cytokines  under  consideration  in  this  study  include  interleukin 
(IL)-1 ,  IL-2,  IL-4,  IL-5,  IL-6,  IL-8,  IL-10,  IL-12,  IL-13,  IL-23,  IL-27,  interferon  (IFN)-y,  and  tumor  necrosis 
factor  (TNF)  -a.  Because  of  the  great  number  of  cytokines  utilized  by  the  immune  system,  some  of  these 
cytokines  must  be  grouped  together  into  a  single  node  so  as  not  to  exceed  computational  limitations.  In 
accordance  with  Folcik  et  al.  (201 1 ),  cytokines  were  grouped  into  either  a  monokine  (MK)  or  cytokine  (CK) 
group.  MKs  represent  the  cytokines  released  primarily  by  the  monocytes  (DCs)  and  CKs  represent  the 
cytokines  released  by  the  lymphocytes  (NK  cells,  Th  cells,  and  CTLs).  MK1  and  CK1  are  pro-inflammatory 
cytokine  groupings,  whereas  MK2  and  CK2  are  anti-inflammatory  cytokine  groupings.  See  Table  4  for  the 
cytokine  groupings  used.  It  should  be  noted  that  the  MK2  node  can  also  represent  the  anti-inflammatory 
dendritic  cells,  however,  since  these  dendritic  cells  have  only  one  input  (Infection)  and  only  one  output 
(MK2),  these  2  nodes  can  be  compressed  into  a  single  node  consisting  of  the  MK2  cytokines. 

Although  most  of  the  immune  cells  exhibit  numerous  levels  of  activation,  since  they  will  eventually 
convert  to  a  state  where  they  release  the  cytokines  of  interest,  they  can  be  condensed  into  a  single  node, 
capable  of  being  downregulated  (-1),  nominal  (0),  or  upregulated  (1 ).  For  example,  Thl  cells  start  as  naive 
germ  cells,  which  differentiate  to  various  levels  of  effector  cells.  Effector  cells  release  cytokines  or  mature 
into  memory  cells.  Since  we  are  only  concerned  with  the  cytokines  these  cells  release,  a  single  node  can 
represent  all  activation  levels  of  Thl  cells. 

Even  though  Folcik  et  al.  (2011)  includes  macrophages,  Tregs,  B  cells,  and  Thl  7  cells  in  their 
model,  these  cells  were  not  essential  in  attaining  any  additional  or  interesting  stable  attractor  states.  For 
example,  B  cells  primarily  release  types  of  immunoglobulins  (antibodies),  which  are  not  measured  in  the 
patient  cohorts  for  GWI  and  CFS.  Therefore,  including  B  cells  in  this  model  would  not  enhance  the  potential 
alignment  with  the  clinical  data.  Also,  some  of  these  cell  populations  included  by  Folcik  et  al.  (2011)  do  not 
have  a  substantial  feedback  role  with  the  majority  of  the  immune  response.  For  example,  Th17  cells  release 
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predominantly  IL-17,  which  does  not  exert  significant  feedback  on  the  other  immune  cells  included  in  this 
model. 


Next,  the  hormonal  inputs  of  cortisol  and  testosterone  were  added  to  the  wiring  diagram  [Fig  2]. 
Cortisol  acts  to  suppress  NK  cell  activity  (Bush  et  al.,  2012),  Thl  activity  (Liberman  et  al.,  2009),  and 
inflammatory  DC  activity  (Zen  et  al.,  2011).  Cortisol  also  has  an  effect  on  Th2  cells,  but  in  a  more  indirect 
way  than  the  other  cells  previously  mentioned.  Cortisol  acts  on  naive  T-cells  to  repress  both  the  Thl 
promoting  transcription  factor,  T-bet,  and  the  Th2  promoting  transcription  factor,  GATA-3,  although  the 
suppressive  effect  on  T-bet  is  much  stronger  than  the  suppressive  effect  on  GATA-3.  Therefore,  if  T-cells 
are  already  committed  to  differentiating,  a  Th2  favored  shift  in  T-cell  development  is  observed  with  long¬ 
term  exposure  to  cortisol,  which  can  be  characteristic  in  chronic  illnesses  (Liberman  et  al.,  2009). 
Consequently,  cortisol  is  assumed  to  have  a  promoting  effect  on  Th2  cells. 

It  has  also  been  demonstrated  that  the  pro-inflammatory  cytokines,  predominantly  IL-1 ,  IL6,  and 
TNF-a,  activate  the  HPA  axis  at  the  hypothalamic  level  by  inducing  CRH  secretion  consequently  leading  to 
increased  release  of  ACTH  and  eventually  cortisol  (Cutolo,  2004). 

In  regards  to  the  male  HPG  axis  effects  on  the  immune  system,  studies  have  shown  that 
testosterone  and  other  androgens  exert  a  suppressive  effect  on  the  IL-6  gene,  thereby  decreasing  its 
synthesis  (Cutolo,  2004).  Testosterone  also  acts  to  enhance  the  Thl  response  and  activate  CTLs 
(Gonzalez  et  al.,  2010).  Cytokines  can  also  exert  negative  feedback  on  the  male  HPG  axis.  It  has  been 
shown  that  receptors  for  the  pro-inflammatory  cytokines,  IFN-y  and  TNF-a,  exist  on  male  leydig  cells  and 
act  to  decrease  the  production  of  testosterone.  As  well,  TNF-a  has  been  shown  to  decrease  the  release  of 
GnRH  in  the  hypothalamus  and  LH  in  the  pituitary  gland  eventually  causing  a  decrease  in  testosterone 
levels  (Foster  et  al.,  2003). 

Not  only  is  there  cross  talk  between  hormonal  axes  and  the  immune  system,  but  also  between  the 


hormonal  axes.  For  example,  the  HPA  and  the  HPG  axes  have  a  negative  effect  on  each  other.  It  has  been 
observed  that  testosterone,  via  its  androgen  receptor,  will  inhibit  the  higher  levels  of  the  HPA  axis  (Viau, 
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2002).  Furthermore,  cortisol  has  been  shown  to  suppress  the  HPG  axis,  inhibit  the  effects  of  testosterone, 
and  downregulate  androgen  receptor  expression  (Mehta  and  Josephs,  2010). 

From  this  more  completely  represented  diagram,  we  can  write  out  the  image  vectors,  ,  for 

each  node  [Table  5],  An  in-house  Python  script  was  used  to  analyze  this  state  network. 

Discrete  Encoding  of  Clinical  Data 

Once  the  simulation  returns  the  stable  attractor  states,  the  predicted  states  are  compared  against 
clinical  data  for  GWI  and  CFS.  Using  a  cohort  of  27  GWI  affected  males,  17  CFS  affected  males,  and  29 
healthy  control  males,  cytokine  levels  were  measured  at  rest  in  IL-12,  IL-la,  IL-lb,  IL-8,  IL-10,  IL-4,  IL-6,  IL- 
23,  IFN-y,  IL-2,  TNF-a,  IL-5,  and  IL-13.  Levels  of  cortisol  and  testosterone  were  also  measured  in  these 
cohorts.  Note  that  only  cytokine  levels  are  being  compared  because  data  on  cell  populations  was  not 
available  from  the  cohort  data.  Also,  the  node  MK27  was  not  included  in  this  comparison  because  levels  of 
IL-27  were  not  included  in  this  study.  The  log  base  2  transform  of  this  data  was  taken  and  the  average  value 
of  these  log  transforms  was  recorded  for  each  cohort.  A  one  tailed  t-test  was  performed  to  ascertain 
whether  there  is  a  significant  increase  or  decrease  in  cytokine  or  hormone  levels  between  healthy  controls 
and  GWI  or  CFS  men.  Ternary  values  are  assigned  to  each  cytokine  based  on  their  p-values.  If  the  relative 
levels  of  cytokines  in  GWI  or  CFS  patients  are  found  to  be  insignificantly  different  from  those  recorded  in  the 
healthy  control  (p-value>0.10),  a  ternary  value  of  0  is  assigned  in  order  to  match  the  model’s  ternary 
system.  If  the  relative  levels  of  cytokines  in  GWI  or  CFS  patients  are  found  to  be  significantly  different  from 
those  recorded  in  the  healthy  control  (p-value  <0.10),  a  ternary  value  of -1  or  1  is  assigned  based  on  the 
relative  difference  between  the  log  transformed  means  of  GWI  or  CFS  patients’  cytokines/hormones  and 
healthy  controls’  cytokines/hormones.  For  instance,  if  the  log  transformed  mean  of  the  GWI/CFS  patients’ 
cytokines/hormones  is  higher  than  the  log  transformed  mean  of  the  healthy  controls’  cytokines/hormones 
and  significance  is  found  (p-value  <0.10),  then  a  value  of  1  is  assigned  to  that  node.  As  well,  if  the  log 
transformed  mean  of  the  GWI/CFS  patients’  cytokines/hormones  is  lower  than  the  log  transformed  mean  of 
the  healthy  controls’  cytokines/hormones  and  significance  is  found  (p-value  <0.10),  then  a  value  of  -1  is 
assigned  to  that  node  [this  data  is  summarized  in  table  6  for  GWI  and  in  table  7  for  CFS].  Since  groups  of 


54 


cytokines  are  modeled  instead  of  individual  cytokines  (ie  MK1 ,  MK2,  CK1 ,  and  CK2),  the  ternary  values 
represented  by  the  individual  cytokines  must  be  aggregated  into  one  ternary  value  representing  the  groups 
of  cytokines.  Taking  the  average  of  the  ternary  values  for  each  cytokine  or  hormone  within  the  grouped 
nodes  and  then  rounding  that  value  to  the  nearest  integer  represents  the  overall  ternary  value  for  that  node 
in  the  model  [Table  8],  These  numbers  are  used  for  comparison  between  clinical  data  and  the  predicted 
values  of  the  model. 
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Figure  Captions 


Figure  1:  Wiring  Diagram  of  the  interactions  between  innate  cells  (yellow),  adaptive  cells  (orange),  and  the 
cytokines  they  produce  (blue).  The  Infection  node  (grey)  is  a  necessary  component  of  the  program,  but  was 
kept  consistently  nominal  so  it  would  have  no  effect  on  the  overall  network  as  per  the  logical  function.  Green 
connections  indicate  an  activating  influence  and  red  connections  indicate  an  inhibiting  influence. 

Figure  2:  Wiring  Diagram  for  the  immune  network  with  hormonal  inputs  from  cortisol  and  testosterone 
(purple).  The  nodes,  Stress,  Infection,  and  Male  HPG  axis  (grey)  were  necessary  components  of  the 
program  required  as  inputs,  but  they  were  kept  consistently  nominal  so  they  would  have  no  effect  on  the 
overall  network  as  per  the  logical  functions.  Green  connections  indicate  an  activating  influence  and  red 
connections  indicate  an  inhibiting  influence. 

Figure  3:  Steady  state  1  (SSI)  as  predicted  by  the  simulation.  Baseline  refers  to  a  value  of  0,  above 
baseline  refers  to  a  value  of  1 ,  and  below  baseline  refers  to  a  value  of  -1 . 

Figure  4:  Steady  State  2  (SS2)  as  predicted  by  the  simulation.  Baseline  refers  to  a  value  of  0,  above 
baseline  refers  to  a  value  of  1 ,  and  below  baseline  refers  to  a  value  of  -1 . 

Figure  5:  Discretized  clinical  data  for  GWI  compared  to  the  two  predicted  diseased  attractor  states  (SSI 
[black]  and  SS2  [grey]). 

Figure  6:  Discretized  clinical  data  for  CFS  compared  to  the  two  predicted  diseased  attractor  states  (SSI 
[black]  and  SS2  [grey]). 

Figure  7:  Sammon  projection  of  the  hamming  distances  between  all  five  states.  Green  dots  represent  the 
predicted  model  attractor  states  (SSO(HC),  SSI,  and  SS2).  Red  dots  represent  the  aggregated  clinical  data 
for  GWI  and  CFS.  Axes  represent  the  aggregated  ternary  (-1 ,  0,  1 )  scores  1  and  2  defined  such  that  the 
observed  Cartesian  distances  between  the  points  in  2  dimensional  space  is  the  same  as  the  hamming 
distances  between  the  points  in  5  dimensional  space. 


60 


Table  1:  Ternary  HIGH/LOW  PASS  operator 


A  B 

B 

-1 

0 

1 

A 

-1 

0 

0 

-1 

0 

0 

0 

-1 

1 

1 

1 

0 

Table  2:  Ternary  OR  operator 


A  B 

B 

-1 

0 

1 

A 

-1 

-1 

0 

1 

0 

0 

0 

1 

1 

1 

1 

1 

Table  3:  Ternary  NOT  operator 


A 

A 

-1 

1 

0 

0 

1 

-1 

Table  5:  Image  vector  equations  for  each  node  in  the  system  shown  in  figure  2  including  references. 


Node 

Image  Vector 

Reference 

MK1 

MK(t  +  1)  =  (DC1(0) 

Folcik  ef  al.  (2011) 

MK2  (DC2) 

MK2(t  + 1)  =  (77*2(0) 

Sia  (2005) 

MK6 

MK6(t  +  1)  =  ( DCl(t))()(TEST(t )) 

Folcik  ef  al.  (201 1 ),  Cutolo 
(2004) 

MK23 

MK23(t  +  1)  =  ( DC\(t )) 

Folcik  ef  a/.  (2011) 

MK27 

MK21(t  +  1)  =  (DC1(0) 

Folcik  ef  al.  (2011) 

CK1 

CKl(t  + 1)  =  (NKs(t)  v  CTLs(t )  v  Thl(t )) 

Folcik  ef  a/.  (2011) 

CK2 

CK2(t  + 1)  =  ( Th2(t )) 

Folcik  ef  a/.  (2011) 

NKs 

NKs(t  +  1)  =  ( MK21(t))()(CORT(t )) 

Folcik  ef  al.  (201 1 ),  Bush 
etal.  (2012) 

CTLs 

CTLs(t  +  1)  =  (MKl(t)  v  CKl(t)  v  TEST(t)) 

Folcik  ef  a/.  (2011), 
Gonzalez  ef  al.  (2010) 

Thl 

Th\(t  +  1)  =  (MKl(t)  v  MK21(t )  v  TEST(t))()(MK2(t )  v 

M^6(r)  v  CK2(t)  v  CORT(t)) 

Folcik  ef  al.  (2011), 
Gonzalez  ef  al.  (2010), 

Sia  (2005),  Diehl  and 

Rincon  (2002),  O'Garra  ef 
al.  (2004),  Liberman  ef  a/. 
(2009) 

Th2 

Th2(t  +  1)  =  (MK2(t)  v  MK6(t)  v  C^2(r)  v  CORT(t))() 

(. MK21(t )  v  CX1(0) 

Folcik  ef  al.  (2011),  Sia 
(2005),  Liberman  ef  al. 
(2009),  Kamiya  ef  al. 

(2011),  Sia  (2005) 

DC1 

DCl(r  +  1)  =  (MK23(t)  v  MK21(t)  v  CKl(t))()(MK2(t )  v 

CORT(t)) 

Folcik  et  al.  (201 1 ),  Zen  ef 
al.  (2011) 

CORT 

CORT(t  + 1)  =  (MK\(t)  v  MK6{t)  v 

Cutolo  (2004),  Viau 
(2002) 

TEST 

+ 1)  =  v  cort(i)) 

Foster  ef  al.  (2003), 

Mehta  and  Josephs 
(2010) 
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Table  6:  Statistical  significance  of  the  one-tailed  t-test  performed  between  healthy  control  (HC)  cytokines 
and  hormones  versus  GWI  cytokines  and  hormones.  Assigned  ternary  values  are  also  listed. 


Cytokine/ 

Average  HC 

Average  GWI 

p-value 

Assigned 

Hormone 

Log  Transform 

Log  Transform 

Ternary  Value 

IL-12 

2.15 

1.68 

0.22 

0 

IL-la 

2.94 

2.21 

0.053 

-1 

IL-lb 

4.06 

3.21 

0.035 

-1 

IL-8 

2.46 

3.85 

0.010 

1 

IL-10 

3.14 

3.00 

0.36 

0 

IL-4 

1.71 

1.57 

0.38 

0 

IL-6 

2.26 

3.10 

0.033 

1 

IL-23 

6.70 

7.69 

0.11 

0 

IFN-y 

1.68 

3.22 

0.022 

1 

IL-2 

2.13 

2.92 

0.10 

1 

TNF-a 

2.95 

3.12 

0.38 

0 

IL-5 

2.22 

2.69 

0.20 

0 

IL-13 

1.40 

2.61 

0.0022 

1 

Cortisol 

-3.09 

-2.58 

0.094 

1 

Testosterone 

8.81 

8.48 

0.033 

-1 

Free 

Testosterone 

6.07 

5.92 

0.17 

0 
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Table  7:  Statistical  significance  of  the  one-tailed  t-test  performed  between  healthy  control  (HC)  cytokines 
and  hormones  versus  CFS  cytokines  and  hormones.  Assigned  ternary  values  are  also  listed. 


Cytokine/ 

Hormone 

Average  HC 

Log  Transform 

Average  CFS 
Log  Transform 

p-value 

Assigned 
Ternary  Value 

IL-12 

2.15 

1.96 

0.41 

0 

IL-la 

2.94 

2.25 

0.19 

0 

IL-lb 

4.06 

4.01 

0.46 

0 

IL-8 

2.46 

3.62 

0.067 

1 

IL-10 

3.14 

3.60 

0.14 

0 

IL-4 

1.71 

2.42 

0.26 

0 

IL-6 

2.26 

2.12 

0.40 

0 

IL-23 

6.70 

8.37 

0.029 

1 

IFN-y 

1.68 

2.69 

0.22 

0 

IL-2 

2.13 

4.26 

0.0026 

1 

TNF-a 

2.95 

3.49 

0.23 

0 

IL-5 

2.22 

1.59 

0.18 

0 

IL-13 

1.40 

1.90 

0.13 

0 

Cortisol 

-3.09 

-1.75 

0.0086 

1 

Testosterone 

8.81 

8.52 

0.074 

-1 

Free 

Testosterone 

6.07 

5.12 

0.033 

-1 
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Figure  1 
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Figure  2 
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Figure  5 
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Figure  6 
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Figure  7 


Sammon  Projection  of  Ternary  Hamming  Distances  Between  States 
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